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BY 
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1. Introduction. The random determinants considered here are those whose elements 
are independent random variables having a common distribution, which is symmetrical 
about zero and (for convenience only) has unit variance. 

They arise naturally in considering the problem of solving large systems of linear 
equations on an automatic computing machine, and have a bearing on the relation of 
the machine precision and the order of the largest system which should be tried. 

With no further restriction on the common distribution than has been given above, 


Fortet [1] has determined the second moment to be 
E(D?) = n! (1.1) 


with D, the nth order determinant in question and £, as usual, the expectation. In 


§2. the fourth moment is determined for the same general conditions. Also, the sixth 
moment for which no general formula has been found is given for n = 1, 2, 3, 4, and 
the simpler cases for two kinds of discrete distribution are displayed. 

The remainder of the paper is devoted to the normal (Laplacian or Gaussian) dis- 


tribution. Here the moments are unexpectedly simple; indeed 


9)! 9 — 9)! 
E(D;") = n! ot? ... e+ 3- 2 (1.2) 


2! (2k — 2)! 


As was first brought to our notice by a reader of an earlier draft, and as has since been 
stated by Forsythe and Tukey [2] to whom thanks are due for the opportunity of seeing 
a draft of their paper, this is a consequence of results of 8. 8. Wilks [4] on the generalized 
variance. Nevertheless we give our original derivation, which takes only slightly more 
space, because it leads more easily to the expression of the probability density as an 
integral of Mellin-Barnes type, from which both series and asymptotic expansions are 
obtained, and because it seems more natural in this setting. 

By symmetry, the odd moments of D, are all zero and the probability density, 
p,(D,,), is an even function. The first four probability density functions are, with z = D, 
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and for x > 0, 


pi(x) = (2r)~'”” exp (—2’/2), 


/ o~-1,-s 
p(x) = 2 'e," 
: Sere r x (1.3) 
p3(x) = (2n)7'”” | r exp (-£ — ‘) dr, 
a 
p(x) = 9 A (2x°"" 


K, is a Bessel function of the second kind for imaginary argument. The result for p.(x) 
is well known. 

2. Moments—general distribution. 
D,, may be written as 


If a,;; is the element of D, in row 7, column J, 


Dp. = 2, he Oyj tis. 78° An;, (2.1) 
with (j; , jo, °** , jn) a permutation of (1, 2, --- n), the sum extending over all possible 
(n!) permutations, and the sign according to the parity (evenness or oddness) of the 
permutation. 

Then 
(2.2) 


E(D:) = E|> + a,;, -** Gn;,]’, ( 


0, that the only contributions are from terms like 


| 
| 


and it is apparent, because E(a;;) = 


Qi1 *** Q, so that 
E(D,;) = n! 
since by convention H(a;;) = 1. 
For the fourth moment by the same reasoning, contributions are only from terms in 
D,, in which elements appear either two or four times. In the latter case, the contribution 


e . ° ‘ ° 4 rer ° * 
of the element to a typical product in D; is ms = E(a;;). The corresponding combina- 
torial problem is that of determining the sum of the weights of all 4 X n rectangles, 

of (2.1)). 


each row of which is a permutation of numbers 1 to n (the numbers 7, , j2, --- . 
The weight of a rectangle is the product of the weights of its n columns. The weight of a 
column is m, if the four numbers in it are all alike. The weight is unity if the four numbers 
are two pairs of like numbers, the numbers in one pair being different from those in 
the other. The weight is zero for any other combination of four numbers. 

Putting the first row in natural order (1, 2, --- n) we may write 

E(D}) = nly, . (2.3) 
Here, y, is a function of m, as well as n, but this need not be indicated explicitly. y, is 
the total weight of all of the rectangles which have the first row in natural order. 

Of several ways we have found to determine y, , the following seems to be the simplest. 
Consider only the rectangles of non-zero weight and classify them according to the 
condition of column n. Either four n’s are in column n, or just two are there (the other 
two in some other column). In the first case, the contribution to y, is m,y, In the 
second case there are three ways of choosing the row for the second n, n — 1 ways for 


choosing the column for the remaining two n’s and the contribution may be written 


‘~~ 


1954) THE DISTRIBUTION OF RANDOM DETERMINANTS 99 


3(n — 1)z, , 2, being the total weight of the rectangles which satisfy any one of these 


choices. Hence 
Yn = MYn-1 + 3(n — 1)z, . (2.4) 


the rectangles whose (n — 1)th and nth columns are n — 1, n — 1, n, 
n and n, n, x, x respectively, where x may be any of the numbers 1 to n — 1. If zis 
n — 1 the contribution is y,-2 . For each of the remaining n — 2 choices of z, it is z,_, ; 
for element » may be removed and the last two columns coalesced to the single column 
x thus repeating the form given above. Hence 


Take for z 


n—Il1,n— 1,2, 
Zn = Yn-2 + (mn — 2)%,-1 - (2.5) 
Elimination of z, from (2.4) and (2.5) shows that 
Yn = (Me +n — Lyn + (mn — 1)(38 — mMy)Yn-2 . (2.6) 
With the boundary conditions, y, = my, ye = mi + 3, this completely determines y, . 
The next two values are: 
y, = m; + 9m, + 6, ys = me + 18m, + 24m, + 45. 


Introducing the exponential generating function 


co 


Y(é) = po y,t'/n!, (2.7) 
0 
it is found by straightforward calculation that, if yo = 1, 
Y(t) = (1 — 0 exp t(m, — 3). (2.8) 
It is interesting to notice that for m, = 0, this is the exponential generating function 


for the number of 4 X n rectangles in which each column has exactly two pairs and is 
also the cube of the corresponding generating function for ‘“‘rencontres” or sub-factorial 
numbers. There is also a purely combinatorial interest in the instance m, = 1. Equation 
(2.8) leads to the explicit expression for the fourth moment 


a + 9 
K(D*) = nly 2 (n k + Im — k + 2) (m, — 3)*. 


k! 


» sixth moment, we have been able to find only the following special values: 


ad k=0 


E(D; > 2'(m, + 15m), 
E(D8) = 3'(m3 + 45m,m? + 30m + 270mi — 90), 
E(D’) = 4![m,(m3 + 90m,mi + 120m> + 1080mz — 360) 
+ 45(17m} + 56m + 144mi — 72m, + 96)]. 

These results are verified by the particular distributions given below in Tables I and 
II. In the first, the elements of the determinant take the values —1 and +1, each with 
probability 1/2, so both the variance and m, are unity. In the second, the elements 
V 2, 0, and 2 with probabilities 1/4, 1/2, and 1/4 respectively, and 


take values 
hence have variance unity but m, = 2. The tables give values of k,p,(x) where p,(x) 
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TABLE I 
Density Functions for Determinants of Elements +1 




















n kn knpn(x) 
z=0 x= 21 z=2" xz = 3.2"1 
1 2 0 ] 
2 2? 2 l 
3 24 10 3 
4 29 338 84 3 
5 214 10744 2505 300 15 
TABLE II 
Density Functions for Determinants of Elements 0, + 2 

nr kn knDy (Zz 

r=QO0 xe = Qri2 z= 2.28/23 L 3.2”/2 = 4,2n/2 
1 2? 2 ] 
2 26 38 12 l 
3 y 1214 327 78 9 3 








is the probability that a determinant of order n has the value xz, and k,, is a power of 
two chosen for convenience of tabulation. Both distributions of course are discrete, 
the values of x being multiples of 2""* and 2””’, respectively, in the two cases, and sym- 
metrical: p(—2z) = p(x). In both cases, as n increases, the density tends to peak at the 
origin. 

3. The average value of an arbitrary function of D, , elements normally distributed. 
Let F(z) be an arbitrarily assigned function of a real variable z such that the various 


averages defined below exist. Later on we shall take F(z) to be one or the other of 
F(z) = |z|*"', R(s) > 1, (3.1) 
ra)= 1 for %<2< 2+ Ae and 0 for other values of z, (3.2) 


where the vertical bars denote ‘‘absolute value of’, and x, x + Az are two fixed real 
numbers. Let ¢,(u) denote the average value (expectation) of F(uD,), so that 


g,(u) = av., [av.. F(uD,)], (3.3) 
where av., denotes an average taken over the elements a,, , «-- , a, of the first row in 
D, , and av., an average taken over the remaining elements a2; , --+ , @,, With a,; , -*-, 


a,, held constant. 
This two-stage averaging is for the purpose of building a recurrence, as will appear. 


Consider a,; , --: , a, to be the components of a vector in n-dimensional Euclidean 
space and let ¢;'’ --- , ¢,'’ be the components of a unit vector having the same direction. 
Form a determinant 7' by taking the components t; ,7 = 1, 2, --- n, to be the elements 


of the first column, t;”’ the elements of the second column, and so on to ¢;"’ where the 
only restriction on the elements ¢;"’, s = 2, 3, --- n is that they be components of a 
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unit vector perpendicular to the n — 1 remaining unit vectors and that their signs be 
chosen so that the determinant of 7 be unity. Such a determinant is always possible. 
The properties of 7 and the rule for multiplying determinants show that 


bi ° On 
Di= DLT =i - . - | = 5D,,B, (3.4) 
| Dar ™ 
where 
by, — (aj, t+ eee + Gs bi, = 0, k= 2, 3, 792 (3.5) 
and 
| bo. bon 
| 
Pe or 
ao 


The inner average shown in (3.3) is equal to the average of F(wb,,B) taken over 

- d,, With the n elements a,, , --~- @;, held fixed. Since b,, and the t;"’ (s,7 = 1, 
, n) depend only on a,, , «++ , a, , they also are fixed. From (3.4) it is seen that the 
elements of B are linear sums, with fixed coefficients, of terms chosen from a2; , «++ Gan. 


Since the a,,’s are normally distributed about zero, so are the b,,’s for the purposes of 


the present averaging process (and for r, s = 2, --- , n). The orthogonal and normal 
properties of the ¢;"’ lead to 
av. b,Dim _ S280 » (3.6) 


where 6,, = 1 if r = l and 0 if r # l. Consequently, in averaging over da; , +++ Gan, 
the elements of B behave like (n — 1)* independent random variables distributed norm- 
ally about zero with unit standard deviation. Thus B is the same sort of random variable 


as D,, except that its order is (n — 1) instead of n. From the definition of ¢,(u) as an 
average we get 
av.; F(ub,,B) = ¢,-:(ub,,), (3.7) 
where av., denotes an average taken over boo, «+--+, Dyn 
From this and (3.3) it follows that 
$n(U) = av.) Gn-1(Ud;;). (3.8) 


At this stage b,, is no longer regarded as fixed but as a random variable by virtue of its 
definition (3.5) in terms of the random variables a,, , --- , @;, . Indeed, writing v for 


b;, , the probability density of v is 


ee ss 
aT (n/2) exp (—v/2), (3.9) 


and therefore (3.8) becomes 


Q-n/2 a @ 


| Gn-1(uv®)y”?"e-"/? dv (3.10) 


¢x(u) = Tid. 
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, (4) in succession starting 


A 


which, in principle, may be used to compute ¢,(w), 
with 9, (u). 

We now use (3.10) to compute the average value of 
the vertical bars denote absolute value. This calls for the form (3.1) for F(z). Since 


D,, \*~* where R(s) > 1 and 


D, = a; and a,, is distributed normally about zero with unit variance, equation (3.3) 
gives 

y,(u) = av., | ua,, [OO = 2 PT (8/2)u"*, (3.11) 
where we take 0 < u < o. Equation (3.11) gives the 0, 2, 4, --- moments of the normal 
distribution upon setting vu = 1 and s = 1, 3, 5, 


The integrals encountered when (3.10) is used to compute ¢,(u), starting from (3.11), 
are of the gamma function type. The general result is 


Ti s—1] ra ssipe Soe 2)T({s + 1] Q) «-- I'({s ot n— 1] 2) 
4 =¢,(l)=2°°"=> eer . (3.12) 
is Pall r'(1/2)T(2/2) --- T(n/2) (3.1 


This equation may also be obtained by transforming the (s — 1)/2th moment of Wilks’ 
generalized variance. 

As in (3.11), the various moments for the probability distribution of D, may be 
obtained from (3.12) ny setting s = 1, 3, 5, --- Equation (1.2) is obtained by taking 
s = 2k + 1, adjusting the fraction in (3.12) to have 2k gamma functions in both numera- 
tor and denominator, and using the duplication formula for the gamma function. 

Let p,(D,,) denote the probability density of D, and let F(z) have the form shown 
in (3.2). The average value of F(uD,) is then the probability that wD, lies between x 
and x + Ar. When Ax — 0 this and equation (3.3) give 


g,(u) = p,(x/u) Ax/u. (3.13) 


A recurrence relation for p,(D,) is obtained when (3.13) is placed in the recurrence 


relation (3.10) for o,(u 


PAT) = 5°21 (n/2) | ain" Ee ae. (3.14) 


/0 


'? exp (—2°/2), the higher p,(x)’s may be computed 


Starting with p,(r) = (2x) 
in succession. The integrations steadily become more difficult and the work has not been 


carried beyond the case n = 4 stated in the introduction. 
4. An integral for the probability density of D, . Since p,(D,) is an even function of 


. (3.12 gives for R(s — 3 


5 El| D, y= | D, ‘pA D,) dD, . (4.1) 
The expression for p,(x) obtained heuristically from (4.1) by Mellin’s inversion formula is 
%s/2)T(iis + DQ) «6% 4 Is _ 2 
Pr(x z | ri ial vse T(s/2)1 = L|/2) at ' ah. . 2 ds, (4.2) 
2771 . '(1/2)T(2/2) --- T(n/2) 


where x > O andc > O. The integrand in (4.2) is an analytic function of s when R(s) > 0 


and the integral converges absolutely by virtue of the asymptotic behavior of the gamma 
] 


functions. This integral is closely related to those given by 8S. Kullback [3] for the dis- 


tribution of the generalized variance. 
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That (4.2) actually gives p,(x) may be verified by induction from the recurrence 
relation (3.14) for p,(a) (the order of integration may be inverted because both integrals 
converge absolutely). 

A series for p,(x) may be obtained by closing the path of integration in (4.2) on the 
left by a large semicircle. This has actually been done by Kullback in treating the analo- 
gous case for the generalized variance. The complication introduced by the multiple 
poles is illustrated by the following series for n = 3. 


= — /9 
px) = (1 ;? > ay | | oe ¢ ie 2) — 5 wh) — (2k — »|) 


x @ (—2x?/2)* 
-32@,ep17 43) 


=a(a+1)--- (a+ k — 1), and W(x) = I’(a + 1)/T(x# + 1). 


Here (a)yo = 1, (a), — 
When k is an integer 


Wik) = —-C+1+ + ; oo - > v(0) = —C, 


9 9 9 
Wik — 1/2) = —C —2log2+2 ove ee (4.4) 
v(—1/2) = —C — 2 log 2, 
where C = .577 


The first few terms of the power series for n > 4 are 


2 n/a zr x2 1—n/2 
Pa\t) = T(n/2)  (n — 2)12 + (n — 2)'IT[(n — 2)/2 
| - toe 2°"x*) + (—3C — 2 log 2+ 3) + > uo: i=) | +--+. (4.5) 
j=4 


It must be remembered that all this assumes x > 0 and that p,(x) is an even function 
of x. 

Equation (4.5) shows how p,(x) behaves for small values of x. An asymptotic ex- 
pression holding for large values of x may be obtained from (4.2) by the method of 
If we assume for the moment that the appropriate saddle-point s = 85 
0, so that the gamma functions may be replaced by their asymptotic 
expressions, we find that (so), = 2°. The proper root to use turns out to be s, = 2°" — 
2. When we take c = 8» in (4.2), the method of steepest descent gives the fol- 


steepest descent 
is far from s 

(n — 1) = 
lowing asymptotic expression: 


n-1)/2¢ (n?—3n+4)/4, (n—1)(n—2)/2n 
x 2 — exp (—na*’"/2) (4.6) 
(2) ~ —ps : _ “). . 
I n1(1/2)1(2/2) --- T(n/2) 


As a check on (4.6), note that it reduces to the exact forms when n = 1 and 2. As 
further check, we note that it satisfies the differential equation (which may be obtained 
by differentiating (4.2) n times) 


a” n_.2-n . — 
qe p(x) — (—1)"2"""p,(x) = 0, (4.7) 
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in the sense of being the leading term of an asymptotic series solution. In fact (4.7) 


may be used to extend (4.6). For example, if n = 3, we can assume 
p3(z) ~ Ax® exp [—327°/2][1 + ax? + ax“? +--+] 
and find from (4.7) that a, = 5/18 and a. = —35/648. 
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ON SOME EFFECTS OF VELOCITY DISTRIBUTION IN ELECTRON STREAMS* 


BY 
8S. V. YADAVALLI 
University of California, Berkeley, California 


Abstract. Based on the Liouville theorem an integral equation is obtained for the 
solution of electron beam problems having a velocity spread. Assuming a rectangular 
velocity distribution (which is justified later) the integral equation is solved by Laplace 
Transforms to obtain the solution of the problems of small-signal velocity modulation 
in a drifting electron stream, and a drifting electron stream initially possessing full 
shot noise in each velocity class. 

It is shown that one obtains from the above integral equation the same results as 
those given by the Llewellyn equations for the case of a single valued velocity stream. 
The problem of finite but narrow velocity spread in the case of an accelerated electron 
stream is briefly considered. 

1. Introduction. To investigate the consequences of treating the electron stream as 
a plasma one can use either one of the two following methods: 

1) treat the electron stream as being composed of a number of beams whose velocities 
are discrete, or 

2) use the distribution function treatment. 

The first method has been used by many; quite recently Bohm and Gross [1] used 
the same to discuss many characteristics of the plasma. The second method has been 
used by Vlasov [2], Landau [3], and recently by Watkins [4] to discuss the behavior 
of a plasma, the effect of a velocity distribution on noise in electron streams, etc. We 
will use below the distribution function method. 

The plasma is completely described by a distribution function of f(r, v, ¢) such 
that —/f(r, v, t)/e gives the average number of electrons in a small range of position and 
velocity at a time ¢ in the phase space whose coordinates are position r and velocity v. 
The behavior of the distribution function is completely specified by the Boltzmann 
equation 

of 
ot 


dv 


_ (of 
+ viV.f + au Vs ~ (if) " 


where dv/dt is determined by the interparticle forces and any external impressed fields. 
The physical assumptions involved in the above are based on the following [5]. The 
forces acting on a particle in a plasma can be divided into two parts. One is the short- 
range force acting on a particle when a close collision is experienced by it during which 
there is a heavy momentum transfer, and this is taken care of by the term on the right 
hand side in (1). The other part is the long-range coulomb force, due to the other particles 
of the system. In the second type of collision, the momentum transfer is small. At any 
instant, a single electron is engaged in a many body collision, and in (1) the effect of 
the same is taken into account as a smoothed out force. This latter force depends on 
the distribution function itself and gives rise to many of the plasma characteristics. 








*Received June 2, 1953. This research was supported in part by the United States Air Force, Wright 
Air Development Center, Air Research and Development Command under Contract No. W33(038)-ac- 
16649 with the University of California. 








106 S. V. YADAVALLI [Vol. XII, No. 2 
That it is so has been shown by Pines and Bohm [5]. In the type of plasmas encountered 
in vacuum tubes the collision term can be neglected, since the collisions would have 
mainly a minor damping effect, although the same procedure cannot be followed in 
the case of a dense plasma. In (1) the value of dv/dt can be obtained by using Maxwell’s 
equations. Also, in the above equation any currents due to the motion of positive ions 
is neglected. To take it into account one can introduce another distribution function 
g(r, v, t) and write an equation similar to (1). If one assumes that positive ions are 
uniformly smeared out one need not consider g(r, v, ¢). 

Equation (1) and the corresponding Maxwell’s equations form a set of coupled, 
non-linear integro-differential equations which are extremely difficult to solve. 

2. Formulation of the problem for a diode type region for small signals. To simplify 
the mathematics we consider below only one dimensional electron streams. 

Consider a parallel plane diode region across which there is impressed a d-c acceler- 
ating potential and a small a-c potential. The distribution function f satisfies the 
specialized Boltzmann equation neglecting collisions 

fy +055 ~ 5 =o (2) 
where x is the space coordinate normal to the diode planes, v is the x directed velocity, 
n = e/m for an electron and £ the electric field. 

We can write 


where V, is the d-c potential and LZ, is the amplitude of the a-c electric field. 
Let 
f(x, v, t) = fol, v) + filz, ve’*’. (4) 


It is assumed that f,; < f, and LE, « FE, . One can then split Eq. (2) into two parts as 


shown below, one containing the d-c terms only, and the other a-c terms as well. 
Ofo OV Ofo , re 
;— eee aE, (5 
Ox Lilie Ox Ov ) 
— Of, , OV, Ofo , Ofo ; 
) = —— —— = wh, = 6) 
Joh FOS Tae oy oy ( 


dropping the products of a-c terms. 
It is evident from the above that f, is a function of wu only, where 


2 2 ¢ r 
“ue =v —2nV,. 


(7) 
The d-c potential V, is measured with respect to the cathode in the case of temperature 
limited flow and potential minimum in case of space-charge limited flow. Similarly the 
velocity u is measured at the cathode or the potential minimum depending upon the 
fact that the flow is either temperature limited or space-charge limited. 
Changing then the variables from x and v to z and u we have, 





vdv=udut+y a dz, dx = dx, 
avd 89 _ 1a in 
wv udw dx Ozu Ox Ou’ 











tN 


+] 
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We also have from Maxwell’s equations 


OE, 1 [ 
— = —— dv 9 

Ox €9 J-a@ fi ‘ ( ) 
where e, is the dielectric constant of free space. After changing the variables according 
to the above scheme, substitution of the same into (6) yields 


-@ Of: _ oy, fo | 
Jy fit ox nb du” (10) 


Integrating the above we obtain 


fi(z, u) = f,(0, u) exp | — ie [ [v(a’, u)]~" az’| 


(11) 


az 


+ [ nu’ 9fo 7,(2’) exp | - ie A [v(é, x)]~* ae] dx’. 


Ou 
Hence, f;(x, uw) can be determined from the above integral equation (11), knowing 
v(x, u), £:(0, uw), fo(u) and E(x). The only boundary condition needed is the value of 
f,(0, u). We can also obtain an integral equation for the current as follows: 


i,(x) = [ vf (x, v) dv = [ uf (a, u) du (12) 
“0 “0 
if the diode is open-circuited for a-c* the total a-c current density is zero, i.e., 
ee ‘ . 
= [t; — jweok,] = 0 
oz [t, jweoE;] , 


yielding 





E, == (13) 


on 
Hence, we obtain after substitution the following integral equation satisfied by the 


current** 


i(x) = | uf (0, u) exp | — ie | [v(ax’, u)}™ az’ | du 
“0 


“0 


(14) 
z ; a © F. ; z i 
+ po [ i,(x’) dx’ | te exp | -ie fy [v(é, u)]° dé | du. 


0 


This is a non-homogeneous integral equation of the Volterra type whose exact solution 
in a general case can be found only by numerical integration. 


*An examination of the Llewellyn-Peterson equations makes it evident that there are a large number 
of problems where the condition of total a-c current is equal] to zero gives the correct result for the a-c 


quantities in the electron stream. 

**It has been found by the author through a private communication that the integral equations (11) 
and (14) have also been obtained independently by L. R. Walker of the Bell Telephone Laboratories, 
and in a somewhat different form by D. A. Watkins presently with the Electronics Research Laboratory 


of Stanford University. 
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3. Drifting electron stream: Simplified integral equation. For a drifting electron 
stream one can write the following simplified forms of (11) and (14): 


/ 


soe , jar. a Ifo cath 
fi(x, u) = f,(0, w) exp (#2) + | nu hi E(x’) exp | ide = #9] ar’, (15) 


Uu “0 


i,(x) = | uf (0, u) exp ( #2) du 
u 


az as ae af o ° zx — , 
+ —_ | 1,(x’) dx’ | ek exp | ea 2 J du 
Jo Ou u 


JW£o Jo 0 


(16) 


which are to be interpreted as follows. The coordinate x is measured from the entrance 
plane of the drift region, and wu is the velocity at the entrance plane. 

Drifting electron stream with small-signal velocity modulation. We will use below the 
integral equation (16) for the solution of a problem treated by Watkins [4] in a different 
manner. 

Instead of using a Maxwellian distribution of velocities as done elsewhere [4] we 
use below a rectangular type distribution for fo(w); the justification for this is given in 
Appendix I making use of the Tchebycheff inequality. 

For the computation of f,(0, wu) the model chosen is as follows. The stream is allowed 
to pass through two parallel grids across which there appears a voltage V,e’**. It is 
also assumed that no transit time effects occur and the grids are perfectly permeable. 


To the left of the gap, we have 


f(u) = & [Stu — u,) — Stu — u, — w)], (17) 
: ” 


where p, is the d-c charge density at the entrance plane, w the width of the distribution 
function, S(u — u,) equals unity for wu > u, and zero for u < u,, u, being the smallest 
velocity of an electron in the drift region defined by its d-c potential. At the right of 
the gap, we have 


PO fol — (u? + 2nVie'*')'7] — Slu — (u? + w? + 2nVie'*')'7J}. (18) 


iP he Fs ‘= 
bs . Ww 


By expanding the step function in a Taylor series and retaining only first order 


terms, we obtain 


V; Pr nV ipo , 
f.(0,u) & —yn — — Ou — u,) ———— §(u — u, — w), 19 
sia aitit. u Ww )+ (u, + w)w (19) 
where 6 denotes a Dirac delta function having the properties 6(a — z,) = 1 forz = z,, 
and 6(x — z,) = Oforz ¥ 2. 


Now we have to solve the integral equation 
i,(z) = uf (0, u) exp (#2) du 
(16) 
wn “a> . — 
* i saad dal | al exp | He) du. 


Jw€o Jo Jo ou Uu 
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To solve the above by Laplace Transforms, introduce the Laplace transform 7,(s) 
defined by 


i,(s) = | i,(a)e~** dx, real part of s > 0. (20) 


0 
Taking now the transform of (16) and using the rule for the transform of a convolution 


we have 





ya} f° eh@, 0 | Ss co : 
i,(s) = | f s+ dele an | 1 ian de ust dale du| , (21) 
nB+j@ 
L(x) = oe | i,(s)e"* ds. (22) 
<7) Jg-jo 


Here, 8 is so chosen that all the singularities of 7,(x) lie in the region where the real part 
of s < 8B. 
We also have, 


aT. * ; 
<2 a, & {itu — u,) — du — u, — w)}. (23) 
ou w 
Substituting (19) and (23) into (21) and intergrating, we obtain 
jwVipon ( iw)( jw ) wp | . 
7 & == ws Ss - 8 naeee-eed ees migene tte, 4 9 24 
u(1 + w/u,) | + u, + u, tw + u(1 + w/u,) (24) 
where w; = np»/€) . One then has to find the poles of 7,(s). These are 
) E w J {Tw w . 4p sae on 
= — a — oo eee . a er a - (29 
2Lu u(1 + w | * 2 [2 ul + w =| 7 u,(1 + w, =} (25) 
Then, = 4 + Be’**, i.e 


(26) 


(t-sartan) +a Meal T 
u, u.(1 + w/u,) u(1 + w/u,) ; 


Neglecting terms of the order greater than one in w/u, , (which is applicable to a 
narrow distribution), we obtain after some manipulation, 


jnV _ w\fw\).. _ W) wed : ~ jek -*)2] 97 
f pi 22)(){sin | (1 2) Us } exp | jo\ 1 u,/ U,) =) 


Notice that as we let w go to zero, we obtain a result identical with that from the analyses 


II? 


1 (2 


of Hahn [6] and Ramo [7]. 
For the sake of comparison the result of Watkins [4] is quoted below in (28). 


wath dstele-Deb—Lob-Deh 2 
Vet Ss u, tt, 4a/ wp ~~ : 4a/ Uu, — joo\ 2c/ u,_|’ (28) 


where, —/,/u, = po of the present analysis, see Eq. (27), and o = eV,/kT, , k being 
the Boltzmann constant and 7’, the cathode temperature. 
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Equation (28) has been obtained using a Maxwellian distribution and a series method 
for solution. One can see from Eqs. (27), (28) and Appendix I that the current amplitude 
obtained with a rectangular distribution is a little less than the one assuming a Max- 
wellian velocity distribution. 

Drifting electron stream initially possessing full shot noise. We next consider another 
simple example to determine the effect of the thermal velocity spread on noise. Consider 
a drifting electron stream initially possessing full shot noise in each velocity class, or 
electrons having a velocity in a small interval around some particular velocity. For 
this computation one needs to find the convection current produced by shot noise in 
each velocity class, and assuming that the noises due to different velocity classes add 
in a mean square manner we obtain the total mean square noise current. This problem 
has been treated by Watkins [4] and by Pierce [8] in different manners from the following. 
Consider then the 7th velocity class, or electrons having a velocity in a small interval 


around 4; . 
1 


f,(0, u;) = [Qew ; po(u ) Af] =u bu — u,), (29) 


i.e., we assume that the input is pure shot noise. The condition on the above expression 


isu, <u; < u, + w, assuming again a rectangular type distribution. We have to solve 


the equation 


; cr? uf.(O, u) P” ar 1 . 
1,(s) = | BANS an | 1 — ple | - — an | (21) 
Jo S++ Jlw/u Jweo Jp OU S + J(w/U) 
(rc ul eu; po(u;) Af]' 2 : ) 
= 4 | ' u—« b(u — u,) due. 
Lo ee gz j\@ u) J 
ih ant q [ poi (uu — u,) — du — u, — wl (30) 
Jwes Jo wis + j(w/u) | J 7 
(Qe Uu; \(u;) Af)! . f Ww 0 lw w -1) =s 
— Rewspolus) Af)" J) 4 2 __ Joo _ (4 d0)(, 4 Jo 
s + j(w/u;) joey u,(U, + w) L\ u,/\ uetw J 


Define now, 
[Qeu;po(u;) Af]'”? = 7; . (31) 


Then we have 


; 2 +e 5 / . a. 
1,(s) jw } Wp ( ia) jw )} } 
= = <le 4. 1 is 4 —><\¢g . kel ee 4 (: 9 
(s L T lu2(1 + w/u,)) | Tt /. oT +. +e J 32) 


m= iy w 
~ | ( ee ee (22) (1 = oY Nt — s,)(s — 8,)(s — s,)]"'. (33) 


2 


At this point to simplify the computations we neglect terms of order greater than one 


in w/u, and write 
-@ Ww Wp w 
wax it (1-4) as%(1-4) 
u, u u, u, 


8 


~~ ... -@ + - Wp 
=? Us J u,’ 
i, (34) 
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Define now 


(35) 


(u; — u,)/u; =e. 
Then the inverse Laplace transformation is taken and after making the indicated approxi- 


mations, one obtains 


i(x) & —i,[exp — jwx/u,]{cos wpr/u,) + jew/wp) sin wpx/u,)}. (36) 


The minus sign appears in the above because of the present notation regarding p) and 
w, . From (36), we find 

. \ 12 “2 2 / 23 2.2 

i(x) |? = tf cos” wpxr/u,) + w/wp)*€ sin” wpr/u,)}. (37) 
Expression (37) is the noise current at x due to an injected noise current in the ith 
velocity class at x = 0. We have 


i; = 2el; Af, (38) 
where J, = p,(u;). u; is the current carried by the ith velocity class. 
We have assumed at the start that the noise due to electrons of different velocity 


classes are independent, i.e., they add in a mean square manner. As is done also by 


Pierce [8] we can write 


i |? > eT, Af{cos? wpx/u,) + (e?)(w/wp)’ sin® (wpx/u,)} (39) 
where, J, = > & I, and 
(?) = (1/I,) DS I; = (1/2) SS Tus — u,)?, (40) 


i.e., (e°) is of the order of the mean square velocity distribution divided by the mean 


Notice that the above result agrees with that of Watkins [4]. 


square velocity. 
The minima of the above expression (39) occur at 
wer/u, = r/2 + ne (41) 

and the minimum value is 


(i?) min Lz {€’)(w/wp)’ el, Af. (42) 


/mnin = 
We next note that Pierce [8] has shown that one obtains the same result when one 
treats the above problem using the Llewellyn-Peterson equations. The assumption in 
the above procedure is that there is no correlation between velocity and current fluc- 


tuations. 
The author [9] has shown recently that the above assumption is valid provided 


one considers that 
a) The fluctuation quantity associated with an electron is independent of the rest 


of the fluctuations, 
b) the fluctuation quantities are identically distributed, and 


c) the process is stationary in time, and observations are made after a steady state 


is reached.* 
4. Solution of integral equation in an accelerating region—single velocity stream. 
It is desired to obtain the Llewellyn form of equations for an accelerating region. To 


*For further details see Ref. 9. 
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do this one has to assume that all electrons at a given plane have the same velocity 
at any instant. We have 


az 


(x) = | uf (0, u) exp | ie | (v(x’, u))~* ax’ | du 


‘ rte F ‘ia 
+— | t(2’) dr’ | = exp E | (v(é, u))™ ae du. (14) 
ou Jat 


JOE 0 /0 


With 


we can write (14) as 


‘ a . . a . oe n Of ? , 
i,(z) = | uf (0, u) exp (—jwr) du + | i,(x’) dx : exp [—jw(r — 7’)] du, 


/0 Jwe ou 

P c+ °s P / / , > ° : 4 
where it is understood that rt = r(x, u) and, r’ = 7’(2’, u). By a proper choice of /,(0, w) 
one can consider either current or velocity modulation or both simultaneously. 


Input current modulation alone. The above integral equation can be written in the form 


t:(2) = 2,(0) exp (—jwr) + | ti(a’)K (ax, x’) dz’, 13) 
where 
. ? os Of» , P 
K(z, x’) = | (2 — exp [—jw(r — 7’)] du (43a) 
Jo \J@€/ OU 
is called the kernel. The solution of (43) can be written as 
i,(x) = 7,(0) exp (—jwr) + | 1,(0) L(x, x’) exp [—jwr’] dx’, (44) 


“0 


where L(x, x’) is called the resolving kernel; it satisfies the following integral equation 


(az, xz’) = K(a, x’) + | Rig... 2iAe’. 2") :de" (45) 


e + "* . ri a . 
as follows from the Fredholm theory [10] of integral equations. Thus, to solve,the integral 
equation (14), one has to solve (45) to obtain the resolving kernel L(z, 2’ 
Now, for the Llewellyn approximation we write 


f,(u) = u "I, du — u,), (46) 


where /, is the d-c current and wu; is the entering velocity, 6(u — u,) is a Dirac delta 
function. 


We then can evaluate K(z, 2’) as follows: 


K(2, 2’) = exp [—ju(r — 7’)] du 


oe. 4 ai re) a a ied ae 
= | fo an exp [—jw(r 7’) IV du. 


JW€o J0 
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Substituting then the value of f, we have 


K(z, 2’) = -|. als | = fexp [—jwr(z, us) + jwr’(x’, u,)]} (47) 


JEU; 


. * . . . . . 
since the delta function vanishes at the two limits. 
Next, one has to find the resolving kernel L(x, x’). At this stage we note that if 


o() = va) + [ o@)K@e, 2’) de’, (48) 
then, 
oa) = va) + | v@)L(a, 2’) ae’, (49) 


L(x, x’) being independent of the form of ¢(x) and (x). We can therefore use the re- 
solving kernel obtained by Knipp [11]: 

L(x, x’) = —{nlo{r(x, us) — 7’(x’, us] exp [—jo(r — 7’)]} [eov(z, u,)o(z’, us) ]'. (50) 
Going back now to the integral equation (43) we have 


i,(x) = i,(0) exp (—jw7) 





™ (51) 
~ | 2,(Q) , exp (—jwr’)} nl o(7 — 7’) {exp [ —jw(7 vated 7’) ]} [eov(x, u(x’, ui} dx’ 
| lo mn Senet % — nl =—_ “ = , Pe -1i o 
0) {exp jor) sale. eas | (r 7’) [v(v’, us] da | 
We note that 
dx’ 

K = * 59 
v(x’, u;) dr’. (52) 

Integrating (51) we obtain 
i(x) = 2,(0)(1 — (nI/2€)(r?/v(z, u;))] exp (—jwr). (53) 


Equation (53) gives precisely the same result as the Llewellyn [12] equations for initial 
current modulation alone. Equation (53) can be written in the form 


i,(x) = E*i,(0) exp (—jwr7) (54) 


as in the notation of Llewellyn. 
Input velocity modulation alone. We have again for the Llewellyn approximation 


fou) = wT, du — u,). 
We assume that the velocity gap is rather narrow. The above equation represents the 
conditions to the left of the gap; to the right of the gap we have 
fot fie’ =u", du — us — v,(Oe’*’), (55) 


where v,(0) is the amplitude of the velocity modulation. To obtain the Llewellyn form 
of equations, we now have to consider that the entrance plane of the electron stream 
is permeable, and v,(0) « u; (the a-c velocity modulation is small in comparison with 


the entrance d-c velocity). 
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Expanding now the delta function in (55) and retaining only the first order terms, 


we obtain 
f(0, u;) & —v,(O0)u"T 6’(u — u,), (56) 


where 6’(u — u;) is the derivative of the delta function with respect to w. 
We use again the same integral equation (43) and the resolving kernel L(x, 2’) 
given by (50). Making use of (56), (43) and (50) as before we obtain 


i,(x) = v,(0)[v(x, u,)]jwrI, exp (—jwr) . (57) 
Equation (57) can also be written as 
. i(x) = v,(0)n ‘jweou; L(x, 0) . (58) 


Equation (57) gives precisely the same result as the Llewellyn [12] equations for initial 


velocity modulation alone, i.e., 
1(2) = F*v,(0) exp (—jw7) (59) 


as in the notation of Llewellyn. E* and F* can be expressed in the forms given in equa- 
tions (53) and (57) after eliminating the space-charge factor ¢ from the Llewellyn co- 
efficients by the use of the equation 
I, : 
n = [u; + v(x) ]2¢[r(a, us]; (60) 
€o 


> 


v(x) = 2n(V.(x), and V,(zx) is the d-c potential at z. 

5. The problem of an accelerated stream with a narrow velocity distribution. The 
successful solution of the simple problems considered in the preceding sections justifies 
the hope that it will be possible to solve more general problems. One can for instance 
make use of the integral equation approach for the case of a finite but narrow distri- 
bution of velocities in an accelerating region as follows. Consider a rectangular velocity 


distribution, 


; Ta a , ' 
Au) = — [S(u —u;) —- SU —- Uz - w)], (61) 
, uw 
where w is the width of the rectangular distribution, S(u — u,) and S(u — u, — w) 
are step functions, 1.e., 
Stu — u;) = 1 for ts 
= 0 for u< U, 


For the case w < u; , one then obtains the following expression for A(z, x’) from 
(43a): 


; . Io w\ Oa - : ; 
K(z, 2’) = ——* (1 -—) — fexp [—jwr(x, ui) + jwr’(x’, u,)]} + A, (62) 
JWEU; u;/ OU; 


where A, a correction term in (62) can be written approximately as 


A= ao (1 _ “) exp [—jwr(x, us) + jwr’(2’, u;)]. (63) 
2 


; 2 
Jwe ou; 


i 
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If we neglect the correction term for the time being, one can see from (62) that the ampli- 
tude of K(x, x’) is lower in this case than in the case of the Llewellyn approximation. 
Finally, there is the general problem considering a Maxwellian distribution of 
velocities instead of the rectangular distribution mentioned above. It appears that this 
problem might be solved by numerical methods employing the above approach. 
Acknowledgment. The writer wishes to thank Professor J. R. Whinnery for his 
suggestions regarding the preparation of the paper and his encouragement during the 


course of the above research. 


APPENDIX 


Let & be a random variable, m its mean value, ¢ its root mean square derivation and 


na number. Then the following called the Tchebycheff inequality* applies to any kind 


of distribution of &, 


P(\t—m|>no) <n”. (A.1) 


The inequality states that the quantity of mass** in the distribution situated outside 
the interval m — no < — < m+ no is at most equal to n-’, and thus gives a good idea 
of the sense in which o may be used as a measure of dispersion or concentration. 

If one assumes that electrons are emitted from the cathode according to a Maxwellian 
distribution of velocities, then the spread in velocity is due to this distribution of 
velocities. Consequently, one can use the above inequality to terminate the distribution 
at some point instead of considering electrons with all possible velocities. The Max- 
wellian distribution of velocities can be written as 


No » (ame do 
nv) dv = (,t)| xp ( st) dv, (A.2) 


where n(v) dv is the number of electrons in the velocity range v and v + dv, ny is the total 
number of electrons, m the mass of an electron, k the Boltzmann constant and T’, the 
cathode temperature. 


We then have 


-) 


| n(v) dv = no 


“0 
(v’) = 2kT./m, (v)? = kT ./m 
ao. = (v’) — wy = (4 — wkT./2m. 
Choosing a value of n = 10, we have no, = 10[(4 — 2)kT./2m]'” we conclude then 
that the number of electrons having a velocity not considered in the above approxi- 
mation is at most one percent. Actually this is a very good approximation which can 
be seen from the solutions of problems in section III of this paper. 

In practical beam-type tubes, electrons are accelerated to a high potential of a few 
hundred to a few thousand volts before they enter the drift region. Consequently a 
little thought would show that one can neglect a small pip of energy ~0.1 volts super- 
posed on an energy of the order of 1000 V. Hence, we can use a rectangular distribution 
whose width is w = no, = 10¢, . 


*See for example S. 8S. Wilks, Mathematical statistics, Princeton University Press, 1943. 
**The probability distribution is here interpreted as a distribution of mass. 
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THE TRANSFER FUNCTION OF NETWORKS WITHOUT MUTUAL REACTANCE* 


By 
AARON FIALKOW (Polytechnic Institute of Brooklyn) 
and IRVING GERST (Control Instrument Co., Brooklyn, N. Y.) 


1. Introduction. In this paper we develop a complete theory of the transfer function 
of general two terminal-pair networks containing resistance, capacitance and self- 
inductance, but no mutual coupling or ideal transformers. These results constitute an 
extension of those obtained in a previous paper [5], for networks containing two kinds 
of elements only. In part, our present techniques depend upon the ideas and methods 
employed in [5] to which reference will be made in the course of the proofs. 

Bott and Duffin [2] have characterized two-terminal networks without mutual 
reactance. The present research is a first step toward solving the corresponding problem 
for two terminal-pair networks. We do not require the results of [2] in this paper. 

tecently several papers [8, 9, 10] have appeared which deal with the transfer func- 
tion of RLC networks. These do not develop a complete theory but are concerned 
primarily with the synthesis of a transfer function up to a constant multiplier. Even 
with this restriction on the multiplicative constant, the synthesis procedures for grounded 
networks described in these papers are all of limited applicability for realizing the general 
transfer function of these networks. Furthermore none of the papers gives the properties 
of the transfer function which are peculiar to general grounded networks and distinguish 
it from the transfer function of general two terminal-pair networks. 

On the other hand the present paper besides characterizing the transfer functions of 
both grounded and general two terminal-pair networks, yields a synthesis realizing any 
multiplicative constant up to the theoretical maximum which is allowable. 

The transfer function A(p) is defined as the ratio of steady state output voltage to 
input voltage in the domain of the complex frequency variable p. It is a real rational 
function which we may write in the form 

A@ =~ KN. xe tert +4 (1.1) 
; D p” + bp” +--+ +5,’ 
where K is a constant and N and D are polynomials which have no common factors. 

We first consider the general grounded two terminal-pair networks (3 external 
terminals) abbreviated 3 T.N. The conditions on A(p) in this case are given in Theorem 
| and are here described in an equivalent form as follows. ‘The poles of A(p) are in the 
left-half plane or on its boundary except that p = 0 and p = = are excluded. A pole 
on the imaginary axis must be simple and have a pure imaginary residue. The zeros of 
A(p) cannot be positive real but are otherwise arbitrary. The range of K is an interval 
0 < K < K, where K, is the minimum value** of the function D(p)/N(p) for0 < p<. 
Conversely when a function A(p) satisfies these conditions a 3 T.N. may be synthesized 
whose transfer function is the given A(p). In the sequel, this synthesis is performed 
assuming an open circuit termination, but the technique may be modified to take account 
of any finite resistive load without any essential change in the above resllts. 

For the general two terminal-pair network (4 external terminals) abbreviated 4 


*Received July 13, 1953. Presented to the American Mathematical Society, Nov. 13, 1951. (Cf. 


Bull. Amer. Math. Soc. 58, 191 (1952)). 


**Ky is a realizable value of K if and only if the minimum value is assumed only at p = 0 or p = © 


or both. 
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.N. the results are given in Theorem 2. They differ from those stated above for a 3 
.N. in that there is no restriction on the zeros of A(p), and the range of K is now 
—K, < K < Ko» where here K, is the minimum value* of | D(p)/N(p) | forO < p <@. 
In the synthesis procedure, it is convenient to distinguish two cases. These are: 
Case I. The poles of A(p) all lie in the interior of the left-hand plane. Case II. A(p) has 
at least one pole on the pure imaginary axis. The synthesis in Case I is relatively simple 
to handle while that in Case IT is considerably more complicated. 

A synthesis procedure in Case II for a 3 T.N., even up to a constant multiplier, has 
not been considered heretofore in the literature. For a 4 T.N., a method due to Kahal 
[8] has appeared which claims to realize the transfer functions of Case II up to a con- 
stant factor as a symmetric lattice. However, both his proof and conclusions are erroneous. 
As a counter-example, we note that the realizable** transfer function A,(p) = 
K(p’ — 0.5 p + 0.5)/(p” + 1)(p + 1) which satisfies his conditions may not be synthe- 
sized as a symmetric latticet for any K # 0. The foregoing example is an illustration 
of a theorem#7 on lattice realization, the proof of which will appear elsewhere. On the 
other hand, Weinberg [9] overlooks Case II completely, and erroneously states [9, 
pp. 37, 53] that every transfer function may be realized as a symmetric lattice. 

2. The grounded two terminal-pair network. 


i i 
ys 


Theorem 1: Necessary and sufficient conditions that a real rational function A(p) given 
by (1.1) be the transfer function of a 3 T.N. are: (i) The zeros of D are anywhere in the 
left-half plane or on the imaginary axis with the origin excluded. 

(ii) At a pure imaginary zero of D, A(p) has a simple pole with pure imaginary residue. 
(iii) The zeros of N can not be positive real but are otherwise arbitrary. 

(iv) m > n. 

(v) The number K satisfies the inequalities 0 < K < Ky where Ky, is the least of the three 
, Lifm = nand of the first two quantities if m > n. If Ky # K, then 


quantities K, , b,,/a, 
K may equal K, . Here K, is the least positive value of « (if it exists) for which the equation 
D — kN = 0 has a positive multiple root. 


Proof: (a) Necessity 

In the case of either a 3 T.N. or 4 T.N. we may write the transfer function (1.1) in 
the equivalent form 

A(p) = Yis/Yu: (2.1) 

where Y,, and Y,, are the short circuit transfer and driving point admittances respec- 
tively of the network, [6, pp. 134-136]. The zeros of D in (1.1) are either the poles of 
Y,2 or the zeros of Y.. which as is well known must lie in the left-half plane or on its 
boundary, [1, Chap. VII]. This establishes (i) except for the exclusion of p = 0. This 
last physically evident condition will drop out of later considerations. 





*Ko is a realizable value of K if and only if the minimum value is assumed only at p = Oor p = 
or both. 

**The transfer function A;(p) with K = 2 is realized by the parallel connection of the networks of 
Figure 5. 

tThis follows from the fact that the fraction (1 — A,)/(1 + A;) has zeros and poles in the right- 
half plane for any K =~ 0, whereas for a symmetric lattice this fraction should be the quotient of the 
two constituent impedances of the lattice. 

ttThe theorem is: Let A(p) be a 4 T. N. realizable transfer function (i.e. one satisfying the con- 
ditions of Theorem 2) belonging to Case II. Define X, = Re {i"-d(A~!)/dp-d"(A-!)/dp"|, n = 2, 3, --- 
Then A(p) can be synthesized up to a multiplicative constant by means of a symmetric lattice if and 
only if at each pure imaginary pole of A the first non-zero value of X, occurs for n even and is negative. 
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The case where D has a zero at p = iw, can only arise if Y,. has a simple zero and 
The remaining possibilities, as is known, 


Y,. has neither a zero nor a pole at p = tw. 
may be excluded by use of the general residue condition 


T1722 — Tho = 0 (2.2) 


for Re(p) > 0, [6, pp 216-218]. Here r,,; , r22 , T12 , Re(p) represent the real parts of the 
short circuit driving point and transfer admittances Y,, , Y22 , Y:2 and of p respectively. 
= iw, of D, A(p) has a simple pole. If Y.. = a(p — tw) 


It follows that at such a zero p = 
+ --. a> 0 then the residue of A(p) at p = tw is evidently Y,. (tw)/a. Since rz. 
(iw) = 0, T1:(two) # ©, (2.2) implies that 7:.(tw.) = 0. Hence Y,.2(iw.)/a is pure imagi- 
nary. This proves (il). 

Condition (iv) is well known and follows from (2.2) and its consequences. To show 
the necessity of the remaining conditions, the 3 T.N. in Fig. 1 is considered upon a nodal 


basis. 














© 

Fic. 1 
Here the ground terminal is taken as node 0, the other input and output terminals as 
nodes 1 and 2 respectively. The remaining nodes are identified so that each branch is 
an R, L and C in parallel. Hence the admittance y,; (¢ # j) of the branch between nodes 
i and j is of the form ap + b + c/p, a > 0,b > 0, c > O. Of course y,;; = y;; , and as 


¥:= Dy; G=1,2,---, 9d, 


i=0 


usual we write 


ji 
where ¢ + 1 is the total number of nodes. Now let y{; = py;; (¢,j = 0,1, °°: ,t,i andj 
not both zero). Then each y/; is a quadratic polynomial in p with non-negative co- 
efficients, and y/; = )>‘-o.;«; yf; (¢ = 1, 2, --- , é). Using the nodal equations of the 


network, it follows* that we may express the transfer function A(p) in the form 


A(p) = A, ‘A, ’ 


where 
, , , | 
Yo Yon ° °° Yor | 
—7 : 1 4 ee —?1 , e rue 
A, Y31 Y33 Yse | _ cop’ + cp Ei ss ge: (2.3) 
—Y¥in ~Yis °° Yir 
; , , | 
es Ee 
—Ybo hast Ube | : :- 
A, =| 432 Yas Ys) = dep’ + dp" + ---+d,,d € 0. (2.4) 
—Yi2 —Yis oh ia Yi 


*Cf. [5, §2(a)] for details. 
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Here A, and A, may have common factors and by (iv) the degree of A, may actually 
be less than s. In [5, Appendix A], it was shown that A, = A, + Aj{ where Aj is the de- 
terminant obtained from A, by replacing the elements of its first column by yi) , — 
Yio ; yio respectively. It was further shown that a determinant of the form A, 
(or A{) is a polynomial in the y{; (¢ ¥ j) with positive coefficients. These two results 


prove that 


O<e, <d,;, (j=0,1,---,8&. (2.5) 

As an immediate consequence of (2.5), it follows that p = O cannot be a pole of 
A(p). Ford, = d,_; = +--+ =d,_, = 0, d,_,-, ¥ 0 in (2.4) now imply c, = ¢,_, = -:: = 
c,-, = 0 in (2.3). This completes the proof of (i). Incidentally, it also follows directly 


from (2.5) that 0 < A(p) < 1 for p in the range 0 < p < ~. (The equality sign may 
hold only if p = 0 or p = &). 

Having established (2.5), the remainder of the proof for the necessity of (v) now 
follows word for word as in [5, §2(a)] starting with Eq. (2.6) there. 
(b) Sufficiency. 

It is useful at this point to introduce a few notations which will simplify the exposi- 
tion. Capital Latin letters (except for A, Y, Z and K) unless otherwise identified will 


always denote polynomials in p with real coefficients. If R = )o%.o a;p' and S = 
p 8;p' then R < S denotes that a; < 8; (¢ = 0,1, --- ,n). If Pis a given polynomial, 


P, and P, will denote the polynomials consisting of the even powers and the odd powers 
of P respectively. More generally, without reference to a given polynomial, the sub- 
script e and 0 will indicate an even or odd polynomial respectively. 

In terms of this notation, the following result is implied by the sufficiency proof of 
[5, §2(b le 
(A) Let S, have only simple pure imaginary zeros and let 0 << R, K S,. Then an LC — 3 
T.N. exists whose transfer function is A(p°) = R./S, and whose Yo = S,/T, with To 
any odd polynomial relatively prime to S, such that S,/T,) is an LC admittance, T, being 
one degree lower* than S, 

For replacing p’ by p in A(p’) gives what we called an R-function in [5, §2(b)]. In 
the RC-network corresponding to this R-function, as given by our synthesis procedure** 
taking p’’’ S.(p'”’)/T.(p'””) as the Y.. , replace each resistance R by an inductance L 
with Z = R and the LC-network required for (A) results. 

Now suppose that A(p) as given by (1.1) satisfies the conditions of Theorem 1. 
We shall construct a 3 T.N. whose transfer function is A(p). The first step in the synthesis 
procedure is to write A(p) in a special form analogous to the special form (#-function) 
used for RC networks. The argument of [5, Appendix B] may be used here to show the 
existence of a Hurwitz polynomialf U (actually the zeros of U may be taken to be 
negative real and distinct except that in the special case where N and D are both even 
functions of p we may use a non-Hurwitz polynomial U having only pure imaginary 
zeros distinct from the zeros of D and from each other) such that in 

A(p) = KUN/UD = G/H 


*This restriction on 7» is easily removed but we do not stop to do this here. 
**Of the synthesis procedures mentioned in [5], the last alternative method in the footnote on p. 120 


which was later given in detail in [4] usually yields a simpler network. 
tIn this paper we shall use the term Hurwitz polynomial in the strict sense, i.e., a real polynomial 
having a positive leading coefficient whose zeros are in the interior of the left-half plane. 
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we have 


0<KG«KH. (2.6) 


For an examination of the proof in [5, Appendix B] shows that in addition to (iii), 
(iv) and (v) of Theorem 1, use was made only of the fact that the zeros of D were not 
positive real or zero, and this is guaranteed by (i) of Theorem 1. This common factor 
technique achieves positive coefficients in the transfer function as is stated in (2.6). 
Recently Weinberg [J. Appl. Phys. 24, 1526 (1953)], citing references to Bode, has 
stated that the use of common factors for similar purposes is well established and almost 
common knowledge. This is hardly justified, for examination of these and other references 
reveals at best a superficial similiarity to our method. Furthermore, in view of his state- 
ment it is surprising that in none of the literature on transfer functions prior to [5] 
including his own thesis [9] has our or a similar method been employed. 

(b,) Synthesis: Case I. 

We now assume that H is a Hurwitz polynomial. Consider H = H, + H, , Hy) = 
pH! . As is well known [7, p 400] H, and H/ are relatively prime polynomials having 
simple pure imaginary zeros and such that H,/H, is an LC-admittance. In view of 
(2.6) we have G = G, + G, , G, = pGi withO <G, KX H, ,0 K Gi « H?. Hence by 
(A) the functions 


A, = Gi/Hi, A, = G./H, 


are LC — 3 T.N. transfer functions. 
We may write 


H 


e Fe pH: ee 
pH: 


8 , 
— ) 73 ——=> 2 = ; 
ptt Ht err en 
where a > 0, a’ > 0,8 > O and F,/H/ and F4/H, are LC impedances, F, and F} being 
of lower degree than H/ and H, respectively. Of course one of @ and a’ is always zero. 
Now according to (A) realize the transfer functions A, and A, by two LC — 3 T.N. 
r, and lr’, whose Y,,’s are Y2,’ = H!/F,, Y2:’ = H./F% respectively. Modify the networks 


lr, and Tr, , without changing their transfer functions, as shown in Fig. 2 to form I} and 














~ 





Cy taking 
- 8 P 
aol tet. Z, = 1+ ’p. 


The Y,,’s of the networks Ij (¢ = 1, 2) are evidently 1/(Z; + 1/Y22’), (¢ = 1, 2) or 


*The method used here is analogous to the last alternative procedure mentioned in [5], footnote on 
page 120 and given in [4]. Similarly, techniques based on the other procedures given in [5] may also be 
employed here. 
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Y,.’s are G)/H and G,/H respectively. Now connect the networks I'{ and I in parallel 
to form a new 3 T.N. T whose Y,, = G/H and whose Y., = 1 so that the transfer func- 
tion A(p) of T is G/H. This completes the synthesis in Case I. 

(b,) Synthesis: Case IT. 

We now assume that H has some pure imaginary zeros. If all the zeros of H are pure 
imaginary, we may apply (A) to obtain an LC-network realization. Otherwise write 
H = P, J where 

P, = []@’ +e) (2.7) 
k=1 
and J is a Hurwitz polynomial. Of course the w, are distinct and different from zero and 
the residue of A(p) at p = + 7w, is pure imaginary. 

The synthesis method here consists of a reduction procedure by which the realization 
of A(p) is made to depend upon the realization of successively simpler transfer functions 
until we finally reach a stage in which we may use the result (A) to effect an actual 
realization. In order to apply this method of reduction, it is necessary that the transfer 
function be written in a particular form //Q (described in Lemma 1) and that a positive 
real function Q/S having special properties (also given in Lemma 1) be associated with 
A(p). In Lemma 2 it is shown how to transform the transfer function into the required 
form as well as how to construct the function Q/S. The essential ideas of the reduction 
algorithm are now described. The justification for the various steps is given in Appendix I. 

Starting with the positive real function Q/S, Q and S of the same degree, we effect 
the decomposition 


Q _ 20’ , Q” 


S S S 


where pQ’/S and Q’’/S are again positive real functions, with Q’ and Q” one degree 
lower than Q. Corresponding to the partition Q@ = pQ’ + Q” we may form the partition 


M = pM’ + M” so that A, = M’/Q and A, = M”/Q” are of the same special form as 
M/Q but of lower degree. 
We have 
S a S’ S S* 
= —7 = bp +7 a> 0, b> 0, 


mpg’ @ Q” 
where Q’/S’ and Q”/S” are positive real functions of the same special form as Q/S. 
Suppose that the transfer functions A, and A, are realized by 3 T.N. I, and I, whose 
short circuit driving-point admittances are Y},’ = Q’/S’ and Y.,.’ = Q’’/S” respectively. 
Then connecting the impedances Z, = a/p and Z, = bp to T, and YL, respectively, as 
in Fig. 2, we get 3 T.N. Tj , P3 whose Y..’s are Y2.” = pQ’/S, Y2.’ = Q’’/S respectively 
and whose transfer functions are A, and A, respectively. Thus the Y,,’s of T'{ and I; 
are Yi.) = pM’/S, Yi. = M”/S respectively. The parallel connection of ['{ and Pj 
gives a network I whose transfer function is 


(pM'/S) + (M/S) __M 





(pQ’/S) + (Q’/S) ~ Q 


and whose Y.. = Q/S. We have thus made the realization of the pair A(p), Q/S depend 


upon the realization of the simpler pairs A, , Q’/S’; A, , Q”/S. 











1954) NETWORKS WITHOUT MUTUAL REACTANCE 123 


As shown in the remarks following Lemma 1, the reduction process may be continued 
until we reach transfer functions belonging to either one of the two following categories 
which we designate as subcases (a) and (8) respectively and whose synthesis we now 


describe. 


(x) Here A = M/Q with Q Hurwitz and of degree 2s. Also0 KM KQ,Q=P.+Q, 
M, = aP, ,0 < a < 1. The required Y,.,. is Q/S with S = P, 
If s > 2, we may expand the LC-impedance P, 7 = P,/pQi as 
P. 
—7 = Bp ar re 
Pp le Qi. 


where B > 0, y > O and F,/Q! is an LC Meth. with F, of lower degree than Q/. 
Now apply (A) to realize the LC transfer function (M,/p)/Q? by means of a network 


lr, whose Y.. = Q//F, . Form the network Ij as in Fig. 2 taking Z, = Bp + y/p so that 
the Y.. of Tj is Q,/P,. . Consider the L-network I, whose series arm is Y, = a and whose 
shunt arm is Y, = 1 — a. The parallel connection of T'{ and Lr, gives a network whose 


transfer function is //Q and whose Y.. = Q/S. 
Ifs = lie. A(p) = [6p + a(p’ + w;)]/[ep + p” + wi] where 0 < 6 < € and the 


required Y,. = 1 + ep/(p* + w;), form the parallel combination of the previous net- 
work I, and an L-network T, whose series and shunt admittances are y,, = 6p/(p> + 
w,) and Yin = (€ — 5)p (p + @) respectively. 


(8) Here A = M/Q withO <« M = M, <Q = P, . The required Y,, is Q/S where 
S = P, + S,, S Hurwitz and of degree 2s. 

Realize the LC transfer function M,/P, -_ (A) to get a network T whose Y,, = 
P./S,. Then form I as in Fig. 2 taking Z, = 1. This is the required network. 

The synthesis method of Case II may es phe be modified to provide an alternate 
synthesis for Case I also. 

A second synthesis method for Case II which is also inductive in nature (but for 
which we give no details) may be described as follows. The transfer function A is written 
as A = rA, + SA, wherer > 0,s > 0,r +8 < 1 and A, and A, are 3 T.N. transfer 
functions. It is easy to show that A is realizable if A, and A, are. We may choose the 
function A, so that it has the same pure imaginary poles as A but its numerator consists 
of even (or odd) powers only. The synthesis of A, is then accomplished easily by 2 
modification of the procedure of Case I. As for A, , it may be taken to have one less pole 
than A on the imaginary axis and the same process is now applied to it until all the 
pure imaginary poles are eliminated, when the results of Case I apply. This synthesis 
results in a network of about the same complexity as the one previously given. 

3. The two-terminal pair network 


Theorem 2. Necessary and sufficient conditions that a real rational function A(p) given by 
(1.1) be the transfer function of a4 T.N. are (i), (ii), (iv) of Theorem 1 and (v’): The number 
K satisfies the inequalities — Ky < K < Ko where Kg is the least of the three quantities 

Kz |, | b/d, |, 1 if m = n, and of the first two quantities if m > n. If Ky # | Ka | then 
K may actually ‘onl + K, . Here K, is that real value of x of smallest absolute ies (af it 
exists) for which the equation D — x N = 0 has a positive multiple root. 


Proof: (a) Necessity 
By the remarks in Section 2 (a) only the proof of (v’) remains. Consider the 4 T.N. 
on a nodal basis taking the input terminals as nodes 0 and 1, the output terminals as 
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nodes 2 and 3, and choosing the remaining nodes as in the 3 T.N. case. Then in the 
notation of Section 2 we have* 


A(p) = (A, — A;)/Az 


with A, and A, as in (2.3), (2.4) and 


Y31 Y32 Y34 “ Y3e 
—Y21 Yoo —Yr “°° Ye 
A; = — Yas —Yi2 Yis ae —Yir ; = cop" + cp’ + id: c, : 
| , / , | 
=n ee). ~Faw ** Ye 


Since A; is of the same form as A, , we have 0 < cj < d; (j = 0, 1, --- , 8) in addition 
to (2.5). As a consequence of these equations, we note that 0 < | A(p)| < 1 for p in the 


range 0 < p <~=, with the possibility that | A(p)| = 1 existing only when p = 0 or 


p =o. The remainder of the proof now follows word for word as in [5, §4(a)] following 


equation (4.2) there 
(b) Sufficiency. 

Let A(p) satisfy the conditions of Theorem 2. We will realize A(p) by showing it 
with A, and A, transfer functions of 3 T.N. If T, and 


can be written as A = A, — A, 
respectively, then the 4 T.N. IT of Fig. 3 


T, are 3 T.N. corresponding to A, and A, 


corresponds to A. 








Fic. 3 


By the argument of [5, §4(b)] together with the remarks in Section 2(b) following 
(2.6), we show that A(p) may be written in the form A(p) = G/H where H +G> 0. 
(b,) Synthesis: Case I. 

If H is Hurwitz and G, , — G, consist respectively of the terms of G having positive 
or negative coefficients, then A = (G,/H) — (G./H) where each of the fractions is a 3 
T.N. transfer function in the realizable form for Case I. 


(b,) Synthesis: Case II. 
If H = P,J, P, =1(p* + w;), J Hurwitz then apply the method** of Appendix II 


*For details cf. [5, §4(a)]. 
**The only modification in Appendix II consists of considering J, + G* instead of G* and J, — G* 


in the argument following (II.1). 
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to get A in the form A(p) = M/P, Q* with Q* Hurwitz, Q* divisible by P, and with 
M, = P, Mt, P,Q%§ + M, > 0, Qt + M* > 0. Denoting the positive terms of M* by 
M! and the negative terms by — M7?’ and using Mj, M%Q’ similarly with respect to M, 
we have A = A, — A, where A, = (Mj + P, M))/P, Q* and A, = (M¢’ + P, M?’)/P, 
Q* are 3 T.N. transfer functions of the form considered in §2(0,). 

We note if one desires the synthesis of A(p) up to a constant multiplier only, a 
procedure based directly upon Case I and having the simplicity inherent in LC synthesis 
techniques is possible. A brief description follows: Write A(p) = KN/D where D = 
P. D’, P, = T(p’ + w,) D’ Hurwitz. Decompose A so that A = (N,/P.) + (N2/D’) 


where V, must be an even polynomial. For sufficiently small positive constants a, 8 


the transfer function A, = a N,/P, may be realized by an LC network T, , and A, = 
BN./D’ may be realized by a Case J RLC network. Then the function [8A,/(a@ + 8)] + 
fa Ao/(a + 8)] = a BA/(a + 8) is realized by first connecting a suitable RLC impedance 


Z, and a suitable LC impedance Z, to T', and I, respectively (similar to the connection 
shown in Fig. 2) and then putting the modified networks in parallel. 

4, Examples. (i) The following example illustrates the theory and synthesis pro- 
cedure for Case I. Consider the synthesis of the 3 T.N. transfer function A(p) = G/H = 
(3p? + p + 2)/(p* + 3p° + 2 p + 3). Calculation shows that K, = K, = 4.468. We 
shall synthesize the given A(p) where K = 3 < K,. Here H is Hurwitz and0 <G < H, 
so that A(p) is in the realizable form of Case I. 

We consider the LC — 3 T.N. transfer functions 





fp?) xx (Go/p) _ __1 2 _ G. _ 3p +2 
AP) = ip) = p+ 2? Aw-7 = 3p° + 3 
Since 
a oe Be ee 
pH. 2 p?2 p +2 
and 
| fy 1 D 
a eo 
H, 3°3 pt+il 
we must realize A,(p*) and A,(p*) with Y,.,’s respectively Y2.’ = 2(p’ + 2)/3p and 
Y33’ = 3(p° + 1)/p. To do this, first realize A,(p) = 1/(p + 2) and A,(p) = (3p + 2)/ 
3 3) as RC-networks whose Y.,.’s are respectively 2(p + 2)/3 and 3(p + 1). This 


(oO P — 
results in two L-networks I, and [, whose series and shunt arm admittances are y,, = 


2/3, vis = 2(p + 1)/3; yor = 3 p + 2, Yoo = 1 respectively. In I, and I, replace each 
resistor R by an inductance L with R = L. We get L-networks ['4, ['% with series and 
shunt arm admittances y* = 2/3p, yi$ = 2(p° + 1)/3p; yX = 3p + 2/p, y% = 1/p. 
Networks ['* and I'% realize A,(p”) and A,(p’) with the required Y.,’s. 

Now connect Z,; = 1 + 3/2pand Z, = 1 + p/3 to 4 and I$ respectively as in Fig. 
2 to form T'{ and I’ whose parallel combination T given in Fig. 4 is the required 3 T.N. 
(ii) The transfer function A(p) = K(p” — 0.5p + 0.5)/(p* + 1) (p + 1), mentioned in 
the introduction as unrealizable by means of a symmetric lattice, belongs to Case II. 
Here K, is found to be 2 and we synthesize A(p) with K = Ky, as a3 T.N. In order to 
apply the procedure for Case II, we work with A(p) in the form A(p) = M/Q = (2p? + 
p’ + 1)/(p? + 1) (p + 1)? together with Q/S = (p’ + 1) (p + 1)?/(p* + Sp* + 5p? + 
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Ve 
1/3 ’ 
= wali zi ay 
5 
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Fic. 4 
5p + 2) which satisfy the conditions of Lemma 1. The synthesis leads us to the network 
formed by the parallel connection of those of Fig. 5 as a realization of A(p). For lack of 
i I £ I 
space, all details must be omitted. 











2/5 Ve 2/s 

——/ PvP. 4} 
2/s 
Ho 

4s 25/g 5/4 

5/4 ' 
Fic. 5 

Appendix I 


Lemma 1. Let M/Q and Q/S be a real rational function and a positive real rational funtion 
respectively satisfying the following conditions: 

(i) The zeros of Q are in the left half-plane including the boundary. 

jii)0o«K<«M<«KQ. 

(iii) Q and S are each of the same degree r + 28, r > 0,8 > 0. 


(iv) Q and S have no common pure imaginary zeros including p = 0. 
(v) 
G P%(iw) (Bw + Bw * + -+: + B,) 
Re S| = SAB Pe) 8 > 0G = 0,1, 5-4) 
cdg a S(t) | 


(vi) Either (a) M, = P.M* ,Q,. = P.Qt , So = P.S&,0OK Mt KQ; or 
(b+) Mp = P.M? ,Q = P.Q3 , S. = PS? ,0K Mi KQ. 
Then polynomials S’, S’’ and a decomposition 
M=pM'+M", Q=7Q'+Q” 


*The polynomial P, is defined by (2.7). 
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exist such that each pair of functions M’/Q’, Q’/S’ and M’’/Q”, Q”/S” satisfy all the 
preceding conditions with r replaced by r—landS = at)’ + ps’ = bpQ” + S” with 
a>0,b> 0. 

Proof: It suffices to consider the case in which vi(a) holds, the other case being treated 
in a similar manner. Our first goal is the decomposition of Q/.S. We note that the zeros 
of S are also in the left half-plane with simple zeros on the boundary because of (i), 
(iv) and the fact that Q/S is a positive real function. Since if 


BO) | _ _ (w") 
haw Ee ~~ | C(iw) |? 
CO) | B(o") 
me Ke ~ | B(iw) |?’ 

t follows from (v) that . 


Re f S| Pi(iea)(Bus” + Bus” * + +++ + 8) | 
Q Q(tw) |" 
Hence any pure imaginary zeros of S must be zeros of P, . Thus we may write S = P? T 
where P! = II'(p’ + w;) of degree 2s’ < 2s, contains all the pure imaginary zeros of S 
and, of course 7' is Hurwitz of degree 2(s — s’) + r. 
[It follows that Q/S may be expanded to yield 


then 


Q = , Oy pe. Pn L, 
S ns px p +w +5 - Pt +7 (I.1) 


where pi ranges over the zeros of P! and a, > 0. Also from ners 


| _ wen pyzasp . Gat” + Ba + >> +8) | 
ti [P..(iw) /P2(iw) | Tia) - 


Let 8! and 87’ (a 1, 2, --- ,r — 1) be any positive numbers such that 8; = 8{ + 67’. 


Let 6, = 8), 8!’ = B, . Form the positive real rational functions pR’/T and R’’/T whose 
numerators are of degree 2(s — s’) + r at most and where* 





[ »R’ '  (Bles?? 2r oe *_e') 

Re &é | = [P.(iw)/Pi(iw)} - le pe ot _—_F) ; (1.2) 
eT pene. Ae te s+: $2) 

Lew | 7 | = [P(iw)/Pi(iw) | - | Tw) P (1.3) 


Then both R’ and R” are of degree 2(s — s’) + r — 1. We have (p R’/T) + (R”/T) 
R/T since the real parts of each member are equal for p = tw and in this case the function 
is uniquely determined by its real part. 
Now split each a in (1.1) so that a, = af + a’, at > 0, af’ > O and write Li/P? = 
o,), Li’/P, = b By (af'p)/(p° + w:). Then 
, 7 
QW 4 e (I.4) 


S S 


> 


La \&P)/ (Pp 


*For the cons struc tion of a positive real function with given denominator and given real part on the 


imaginary axis see, "13, pp. 97-98]. 
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where 
pt)’ = LiT + pR’'P? , Q” = Li’T + PIR”. (1.5) 
Q’ and Q” are each of degree 2s + r — 1. The desired decomposition of Q/S is given by 
(1.4) and (1.5). 
Before proceeding to obtain a corresponding split of @, we must first show that 
P, divides both Qj and Q?’. As a preliminary to this we need the result 


1, =P,U,. (1.6) 
We have 
- 2] —_ Qel ies) So(— tie) + Qu(Gs) So( =i) 


* S(iw) 


_ Pfiw)|Q*(iw)S.(—itw) + Qo(tw) S7(— iw) 

= S(tw) ; 
where we have used vi(a) to get the last form. In view of (v) it follows that P,(7w) (and 
a fortiori P! (iw)) divides [Q* (tw) S, (—iw) + Qo (tw) S% (iw)]. Now P? (iw) divides 


S, (—iw); but because of (iv) P/ (tw) is relatively prime to Q, (iw). Hence S* (iw) is 
divisible by P? (iw). Thus Ty) = S,/P? = P, (S%/P) which is (1.6) with U S3/r.. 

Coming back to Q’ and Q”’, we show now that 
Q:’ = PC,. (1.7) 


In view of (1.3) and using an argument similar to that given above in the proof of (1.6), 
; co 5 


we find that [R2’ (iw) T, (— iw) + Ri’ (iw) To (— iw)] is divisible by [P, (iw) /P! (iw)]’. 
Since T is Hurwitz and by virtue of (1.6) we have P, dividing 7, but relatively prime to 
T. . Hence P,/P! divides R?’. Using this result in the second equation of (1.5) and again 
noting (1.6), we find that (1.7) holds. A similar argument using (1.2) and the first equa- 


tion of (1.5) shows that 
Q=P.D. (1.8) 
We are now ready to decompose M. Equating the odd and even parts of the numerators 


in (1.4) and deleting the common factor P, , we get 
Qo = pQ: + Qi’, QF = pDi + C.. (1.9) 


Since 0 < M, K Q, and0 « M* « Q% by (ii) and (via), we may obtain* polynomials 
Ni, No’, Bi, Bt’ such that N§ + No’ = M,, Bi + BY’ = Mt andO0XKNj XK pQli,o<« 
Ni’ KX Q0',0 KX Bi< pDi, 0K BY’ KC, . Now define M’ = (No + P, Bi)/p, M" = 
Ni’ + P. BY’. Then M = pM’ + M”. Since B is divisible by p’, M’ is actually a poly- 
nomial. Also P, divides both M/ and M2’. 

Finally we obtain polynomials S’ and S’’. In pQ’/S, Q”/S, take out the zero at 
p = 0 and p == respectively in the usual way. Then we may write 


Y ose 


S a S’ S A 
7 = —, w=! Pv d, 
2Q . + Q’’ Q yp + Q”? a> 0, b>0 


*We are here using the following result whose proof is straightforward. If R and S are any poly- 
nomials with real coefficients and0 <«R<«S,0<«S,,0<«S,, S,; + S. = S then there exist R; and 
R, such thatO < R; K S: ‘ 0< R, <K S, ’ R, + R, = Rf. 
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Here S’/Q’ and S’’/Q” are positive real functions. The relations 
S = aQ’ + pS’ = bpQ” + 8” (T.10) 


hold, and together with (1.7), (1.8) and vi(a) imply that P, divides both S{ and S%’. 
Also we have Re (p = iw) [S’/Q’] = Re (p iw) [S/pQ’]; and since by (1.5) and (1.2) 


” Es oie 3 = Pig) - Bae + Blu” + +++ + Br") 
ae s S(iw) |? 


S 7 


it follows that 


, y 2r-2 y 2r-4 ee , 
Re ka = P%(iw) - ln i : 


Ss | Sw) 


' 


Similarly, we get 


QO”) _ py. Bie + Bie + + + 8 
Re E | = Pifie) - Bre + Bee tt Br) 
, A S (tw) 


; — > 


By the last two equations S’ and S” are both of degree 2s + r — 1. We have now obtained 
the required pairs of functions M’/Q’, Q’/S’ and M’/Q”, Q”/S”. The preceding proof 
contains a verification that they satisfy all the conditions of Lemma 1 except (iv), (with 
1). In particular the first pair satisfies vi(b) and the second pair vi(a). 
To establish (iv) for Q’ and S’ note that by (1.10) a common factor p* + a’ of Q’ and S’ 
must divide S and hence P!. But then by the first equation of (1.5) it must divide Lj 
A similar argument proves (iv) for Q’” and S”. This completes the 


a replaced by 


which is impossible 
proof of Lemma 1. 

The reduction procedure of Lemma 1 may be continued for each pair of functions 
until we reach pairs for which r = 0. Let then A = M/Q and Q/S satisfy the conditions 
of Lemma 1. with r = 0. We now investigate what this implies as regards the form of 
M, Q and S. First suppose condition vi(a) holds. We may take Q = P, + Q, . Since 
S is of degree 2s and P, divides S, we must have S, = 0. But then using (v) we find that 
P. (iw) S, is divisible by P? (tw). Thus S = S, = dP, ,d > 0 and without loss of 
generality we may take d = 1. In view of (i), (iv) and (v), Q is Hurwitz. Also by vi(a) 
M, ar. ,0 <a s f. 

Next suppose condition vi(b) holds. Then Q, = 0, and by (ii), 44, = 0. We may take 
S P, + S, , and as in the preceding case show that Q = Q, = P, , and that Sis 
Hurwitz. 


Appendix II 

Lemma 2. Let A(p) = G/H where®*0«&G KH, H = P.J, J Hurwitz, and suppose A(p) 
satisfies the residue condition (ii) of Theorem 1. Then functions M/Q and Q/S can be 
constructed having the properties slated in Lemma 1 and such that A(p) = M/Q. 
Proof:** lf P, does not divide J, then we begin by constructing a Hurwitz polynomial 
F such that [FJ], is divisible by P, . Let J, = P! J* where P{ =. II’ (p” + «;) contains 
the common zeros of P, and J, . 

*The polynomial P, is defined by (2.7). 

**Our construction will yield a pair of functions which satisfies vi(a) of Lemma 1. By a parallel 


procedure we can get a pair which satisfies vi(b). 
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Let iw, be a zero not included in P! and suppose J(iw,) = a + 18, a ¥ 0. Form the 
LC-impedance 7fF(/F! with n constant and Fj and F? relatively prime polynomials 
subject to the following further conditions: (i) P/ divides F; (dia). If 8 = 0, p’ + w. 
divides F!, » = 1; (iib). If 8B ¥ O choose the remaining zeros of Fj and the zeros of F? 
such that sgn [7F! (iw,)/F6 (iw,)] = sgn (8/a). Take n = taF? (tw,/BF 6 (tw,). 

It follows that » > O and F’ = Fi + 7 Fj is Hurwitz. One readily verifies that [F’/], 
is divisible by (p* + P!. The process is continued with F’J to determine an F’”’, ete. 
until all of the w, are exhausted. Then F may be taken as the product* F’ F” -- - 

Now consider A(p) = GF/P.JF. Since P, divides [JF], it follows from the residue 
condition (ii) of Theorem 1, that P, divides [GF], . We may therefore suppose that 
function A(p) of Lemma 2 has already been prepared so that G, = P, Gt, J, = P.J%. 
Of course 

0<G, «P.Jo , 0x G, = FG! << PJ. . (11.1) 


If also 0 «< G* « J, , then we have one of the special forms of the transfer function 


required for Lemma 1. In the contrary case, we will bring this about by means of a 
alate G* = (PJ. = Gy P. : 


suitable common factor. Consider the polypomials G* and J, : 
In view of (II.1), both of these even polynomials have positive leading coefficients and 
have no positive real roots. Hence by [5, Appendix B] (with p replaced by p”) there exists 


a polynomial EF, = II(p’ + 6.) with the 6, real, distinct and different from the + w, and 
zero, such that 
0<£E.G?* <E,J,. (11.2) 


We now choose as our common factor the Hurwitz polynomial 
B=P.P5+ EE; 
where P/ and E? are any polynomials making P,P{/E,.E! an LC impedance with P,P% 
relatively prime to E.E!. 
Now define 
M=GB, Q*=JB, Q=P.Q*, 
so that A(p) = M/Q. Then we have 
M, = P,|G,.P: + GtE.E?), Q? = P.PitJo + PLJE LE . 
Hence P, divides Q* and by (II.1) and (II.2),0 «< M,/P, M* < Q*. Thus M/Q is 


in the form required by Lemma 1. 
We will now construct the function Q/S mentioned in Lemma 2. We note that Q* 


is Hurwitz and let its degree be r. Choose any positive numbers 6; (k = 1, 2, --- , s) and 
6, (k = 1, 2, --+ , r), and form the positive real rational function** 
S : 5ip V 
ne s+ (11.3 
o-~ trtate 
where 
V Baw” + Bw” + +: +8, A(w 
Re | | o es et Hits a a = = . (11.4) 
ete ar Q* (iw) Q* (iw) 


*A more complicated algorithm which handles all of the w’s simultaneously and results in a poly- 





nomial F of much lower degree can also be given. 
**Cf. footnote regarding Eq. (1.2). 
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Then S and Q are both of degree r + 2s. It follows from (II.3) that S and Q have no 
common pure imaginary zeros. Also since P, divides Q* and 


i dip 
S = P,Q* =; + P, 
Q LF + w, + 


we have that S, = P,S%. From (II.3) and (II.4) we get 


Q| _ Pi(iw) Ow) 
~ H = T Su) ? 


p=iw 


Thus Q/S has all the properties required by Lemma 1. 
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NUMERICAL METHODS IN DIGITAL REAL-TIME SIMULATION* 


BY 
H. J. GRAY, JR. 


University of Pennsylvania 


1. Introduction. Under the sponsorship of the Office of Naval Research, Special 
Devices Center, the University of Pennsylvania conducted a study to determine the 
feasibility of building a digital real-time simulator for a system described by thirteen 
first-order non-linear differential equations. To achieve this, it was found necessary to 
improve either the speed of the machine [1] or the quality of the method of solution. 

The present paper is an introductory one and deals with the latter problem and 
especially with the build-up of errors. Essentially the paper falls into two parts. A. 
Finding the characteristic equation associated with the difference equation by which the 
differential equation is replaced for purposes of numerical solution [2] B. Solution of 
the characteristic equation by conformal mapping. A diagram is given showing the 
method of solution in a case of practical interest. 

2. Separation of round-off and truncation errors. This paper, being introductory 
in nature, will consider only linear constant-parameter systems of ordinary differential 
equations. Assume that the differential equations have been replaced by linear difference 
equations. Round-off errors may be treated as perturbations introduced at each step. 
Since the system is linear the propagating effect of a perturbation satisfies the same 
difference equations as the precisely computed solution and, if desired, the propagating 
effects of the individual perturbations may be summed and added to the precisely 
computed solution to obtain the actual computed solution. 

3. Derivation of the characteristic equations. With the knowledge that the precisely 
computed solution and the round-off errors satisfy the same finite difference equation, 
it becomes important to study the properties of these difference equations. 

Consider the system of linear constant-parameter differential equations 


dzx,;/dt = z. Q;;2; 3 $= 1,2, --- ,k. (1) 
7=1 
Define the matrices 

X = column (2, , 22, °** , %); 
¥ 
| Ay, Aye ig ig ay. | 

: 

A= | 
| 
L xy Axe ewe Axr_) 


For simplicity in exposition, A will be assumed to have real elements and distinct 
characteristic values. 
A class of quadrature methods which have proven to be useful in real-time simulation 


*Received May 26, 1953. 
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involve in one time interval, when applied to Eq. (1), the sequence of computations 
given below. (Methods [Ova , Cra]) 


N 


M 
Onn: (ad = Do alae; +h dD, bjo(dz./dt)2,; , 


i=1 7=1 


R Q-1 
Cro: (2)2’ = Doa,(edel; th > b;(dx,/dt).; , 


i=1 i=0 


where 
k 
(dx;/dt)t = Drafe)i; i= 1,2,-+-,k. 
j=1 


A constant time interval, h, is assumed and the ajo , bj. , a;, , and b;, are assumed to be 
known constants as are N, M, R, and Q. The symbol (z;); denotes “‘the ith variable x 
evaluated at time ¢ = jh, j an integer.’”’ The primes and double primes serve only to 
identify the two classes of the numbers (z;); , namely the results of the operations Oy 
and Cro . 

Equations (2) when expressed in matrix form become 


N M 
fini wow 2 eee, ee? OA, 


7=1 7=1 
(3) 
R Q-1 
Cosi EP mw F430, +4 2, t,A4XL,;, 
7=1 7=0 
where X’ , X/’ are column matrices with elements (x;)/ , (7;)}’ respectively. 
Proceeding in standard fashion [3], assume 
X/ = Re X35 exp (nw), Xi’ = Re X34’ exp (nw), (4) 


where X{ and Xj’ are column matrices having complex elements and w is as yet an 
undefined complex number. Substitution of (4) into (3) yields, after elimination of XQ’, 


Q-1 R -} N 
i[ >. b;. exp (—iuy || >> a;. exp (—jw) — | | jo exp (—u) | 
1=0 i=1 


7=1 


M 

+E—h b> bio exp (jw a}, = (”, (5) 
7=1 

where £ is the unit matrix and 0° is the null matrix. 

For a non-trivial solution of Eq. (5) to exist, the determinant of the quantity in the 
curly brackets must vanish. This determinant is found to be the characteristic determi- 
nant of the matrix A, with its characteristic values \ equal to functions of exp w. Simplifi- 
cation of the rather formidable expression obtained in this way yields the characteristic 
equation 


z= f(a + Bf)", (6) 


where z = )h, p’ = exp ( — jw) (see section 4); 


R Q-1 = 
f= E ~ va, || S ba | 


N M 
a = > ajp’, B= >> djop’. 
j=1 


7=1 
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Similarly it can be shown that for repeated closures using Crg , z = £ if | bo.hAA| < 1; 
for Oy used alone, z = [1 — a]8™'; and for the Runge-Kutta method [4], 


4 
expw = >> 2'/jl. (7) 
j=0 


4. Solution of the characteristic equations. Detailed specifications on the quality 
of simulation will not be given here. However, the accuracy requirements are generally 
considerably looser in simulation than in other scientific computations such that the 
conformal mapping method of solution given below has been found to be entirely 


adequate. 
Any general solution x; of Eq. (1) may always be written as a linear combination 
of the functions exp (A,t), 7 = 1, 2, --- , k where the \, are the (distinct) characteristic 


values of the matrix A. The X, are called the “complex frequencies” or sometimes just 
“frequencies” of the system described by Eq. (1). By exploiting the parallel nature of 
the theories of linear constant-coefficient differential and difference equations, it is 
possible to define and talk about the “frequencies in the computed solution’’. Hence, 
determination in the computed solution of the frequency having the largest real part 
enables a determination of at least whether or not a given computed solution of Eq. 
(1) will decay to zero (stable) or grow indefinitely (unstable). This is the source of the 
name ‘‘stability charts’ given to the conformal maps to be derived. 

Equation (7) yields as roots, given an unique value of z, a set of values of w differing 
by integer multiples of 277. However, (see Eq. 4) since exp n (w + 27il) = exp nu, 
n and | being integers, it is sufficient to consider only one value of w. That value satisfying 
the inequalities —z < Jmw < z will be chosen and will still be denoted by the symbol 
w in the following. 

Let 

w = ph. (8) 
Equation (8) defines u, the frequency in the computed solution corresponding to A, the 
frequency in the true solution. 

Equation (7) being in a form in which it can be considered solved will not be discussed 
further. The remaining characteristic equations are written in the form z = f(w). How- 
ever, they may be put into the form 

Aor" + Ay" "+ --- +A, =0, 
A, ~ 0, 


i =i 


p™' = exp ju, (9) 


II 


where the A, are first degree polynomials in z = Ah. Corresponding to a chosen value of 
\, Eq. (9) has n roots r, which when arranged in order of decreasing magnitude are [5] 


Pt a ee ae” = » Tn ° 
Since w = log, r, one has w, , w. , --- , w, Where Re w, > Rew, > --- > Rew, 


and —r < _Imw, < 7m. 
It is evident from the remarks made earlier about stability that w, is of interest. 


A number of methods of finding w, , given A, were investigated including numerical 
extraction of roots and expansion of roots in a power series in \h but these were dis- 
carded in favor of the conformal mapping method for the following reasons. 
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(a) w, as furnished by the conformal maps (stability charts) is sufficiently 
accurate for the simulation application. 

(b) The computation required to prepare a stability chart is relatively routine 
and the amount of it is not much greater than that required for one root extraction 
by standard methods. Furthermore, the stability chart effectively solves all root 
extractions of practical interest for the given quadrature method. 

(c) The stability charts furnish a convenient means of comparing the merits of 

the different quadrature methods. 


In order to simplify the exposition, the following assumptions will be made con- 
cerning the nature of the characteristic equation, z = f(w). 

(a) zis a single-valued analytic function of w except at a finite number of poles. 

(b) w is a multi-valued function of z. The restriction —1r < Imw < aw makes w 


assume only a finite number of values for any one value of z. 


The stability charts are obtained by mapping the w plane into the z plane using the 
characteristic equations. Essentially what this does is to determine one root for each z. 
Whether or not this is the root, w; , is determined as follows: given z = f(w) it is well 

’ b] a) 
known that there may exist for a fixed value of z = z, , a root w, which has the same real 
part as w, . The trio of points z, , w, , w, will be defined to be a “point on a branch 
contour’’. Hence, a “‘branch contour”’’ will be the locus, if it exists, of “‘points on a branch 


contour” in both the z and w planes. Furthermore, every branch point lies on at least 


one branch contour. 

Fix z = 2 and determine the roots w of f(w) = 2; 
distinct). Denote the point in the w plane defined by the w, as w, . 

It can be stated that if a point w*, which lies inside a simply connected region con- 
but not containing any points on a branch contour, is mapped into the z 
then one root of 2* = f(w) is w*, and for all other roots w, , Re 


(assume that all the roots are 


taining w, 
plane by 2* = f(u*) 
WwW, < Re w*. 


This follows directly from the definition of “points on a branch contour’ and the 
theorem [5] that the roots w of z = f(w) are continuous functions of z. 
y v 
Z PLANE W PLANE 


Ws 


Wi 
Cc 








x 





Fic. 1. Determination of points on a branch contour 


To find points on the branch contours, let w = u + i follow C in Fig. 1. Its image 


z = x + zy may follow C, . The crossover point z, and its images on C lie on a branch 


contour by definition. 

The characteristic equation, z = f(w), for [O,, , C,,] with coefficients given in (10) is 
mapped in the rough in Fig. 2. The contours in the w plane are the u = constant lines, 
with v (and u) taken at intervals of 0.1. Symmetry and periodicity restrict the required 
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| 


= ame : 





us~ 





-0.5 -0.4 -03 -02 -01 O O1 02 03 O04 05 x 


Fic. 2. Preliminary mapping for [Ou , Ca] 


range of v to0 < v < x. Forz = 0, it is found that w} = 0. The w plane maps into the 
z plane five times, but since only the dominant root, w, , is of interest, four of the five 
coverings of the z plane have been discarded (some of them in the sense of not even being 


computed) to obtain Fig. 3. 


Ou 2 Qo = —10/3 Cai Qi. = 48/25 
dey = 6 a, = —36/25 
azo = — 2 a; = 16/25 (10) 
Qo = 1/3 a, = — 3/25 
bio = 4 b. = 12/25 


5. Computations. One way of using the charts is as follows. Given the system of 
differential equations, determine the frequencies, A, in the system. Spot the points 
z = dh on the stability chart for the quadrature method to be used and read off the 
values of w = wh, a procedure rigorously correct only for linear constant-parameter 
systems. 

The process will be illustrated by formal application to the linear varying-parameter 
system described by 

dx/dt = y 
(11) 
dy/dt = —3y — (6.25 + 2.4t + 0.362’)x + f(d) 
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0.8 
0.7 
0.6 
0.5 


0.4 







Vv. 
0.3 


Veo 
v=2.9 
V=3.0 


V=3.1 
-0.6 -05-04-03-02-01 O Of O02 O03 O4 O5 X 


Fic. 3. Stability chart for [Os , Ca] 


0.2 


0.1 


whose solution was computed for two cases as follows using method [0,, , C;,] (Eq. 10). 


(a) h = 0.1, f() = + 0.000rxx where the xxx were random numbers, and 
xz; = 1,y; = 15 forj = 0, --- , —3. The results are given by the h = 0.1 points 
of Fig. 4. 

(b) h = 0.02, f(t) = O and starting conditions as in (a) above. The results are 
given by the h = 0.02 curve of Fig. 4 and serve to approximate the true solutions 
of Eq. (11). 

The matrix A is 
0 1 
A= ; 
—(6.25 + 2.4t + 0.362’) —3 
and its characteristic values are 
X= —1.5 + i(2 + 0.60). 
For h = 0.02, 2 = —.03 + 7 (.04 + 0.012¢) which for ¢ < 5 yields nearly the same 


value for w (Fig. 3). However, for h = 0.1, 2 = —0.15 + 7 (0.2 + 0.06) crossing the 
u = 0 line at ¢ = 2. The expected unstable behavior appears shortly thereafter (Fig. 4). 
6. Conclusions. The introduction to a systematic and convenient method is de- 


scribed for determining the stability of a system by operating on the frequencies of the 
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X(t) 





0 1 ~ t 3 4 5 


Fic. 4. Computations with [Ou , Ca] 


system. The procedure is rigorously established for linear constant-parameter systems 
only, even though an example is given for a linear varying-parameter system. 
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APPENDIX 


The stability chart for [O,, , C,,] together with charts for thirty-two other quadrature 
methods was computed directly from Eqs. 6 on the automatic digital computer at the 
University of Pennsylvania Computing Center. Specifically, for [O,,, C,,], one has 


a= —10/3p+6p’ —2p'4+ 1/37", 
B= 1p, 
¢ = [1 — 48/25 p + 36/25 p® — 16/25 p® + 3/25 p‘][12/25]"’, 


z=fla+ pt)"; p= exp(—w). 
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LOWEST EIGENVALUES OF NEARLY CIRCULAR REGIONS* 


BY 
HERBERT B. KELLER AND JOSEPH B. KELLER 


Institute of Mathematical Sciences, New York University 


1. Introduction and formulation. Consider a plane, simply connected, closed, domain 


bounded by the curve 


r= KO), (1) 
where R(@) = R(@ + 2z). It is well known that the problem 
(V7? + dA‘)u(r, 6) = 0, r < R(6); u(r, 6) = 0, r = R(6), (2) 


for sufficiently regular curves (1), has a non-trivial solution if and only if \” has one of an 
infinite sequence of positive real values, the eigenvalues. In many problems of mathe- 
matical physics (e.g. waveguides and vibrating membranes) the smallest such value, )’, 
called the principal eigenvalue, is of importance.t In the present paper a method for 
approximating the principal eigenvalue is considered. 

The principal eigenvalue of any domain lies between the principal eigenvalues of the 
inscribed and circumscribed circles, as is well known. If the principal eigenvalue of the 
circle of equal area is used instead of that of the circumscribed circle, the same statement 
applies, as is also well known. The latter statement is sharper in that the equal area circle 
furnishes a better (i.e., greater) lower bound than the circumscribed circle. 

In this paper we attempt to improve the lower bound, that is, to get a still greater 
lower bound. We do this by introducing a family of functionals of the domain which 
depend monotonically upon a parameter a. For a = — ©, 2, +© this functional yields 
the principal eigenvalues of the inscribed, equal area, and circumscribed circles respec- 
tively. We then conjecture that there exist values of a which furnish better lower bounds 
than does the equal circle. For almost circular boundaries we find that a = —1.779 --- isthe 
best such value. 

It seems probable that when a = —2, in which case the functional is quite simple, a 
good approximation which is not necessarily a lower bound, is obtained for \*. As an 
example of the utility of the result, the exact eigenvalues of a rectangle with arbitrary 
ratio of sides is compared with the approximate expression, given by an integral, for 


various values of a. The agreement is within 3% for all rectangles when a = —1 and 
even better for nearly square rectangles when a = —2. 


2. The a-mean radius.** The boundary (1) is given parametrically in terms of 
rectangular coordinates with origin at the pole and z-axis on the ray 6 = 0 by 


x(0) = R(@) cos 8, y(@) = R(@) sin 6; = 8 <= De. (3) 


*Received July 11, 1953. This work was performed at Washington Square College of Arts and Sci- 
ences, New York University, and was supported in part by Contract No. AF-19(122)-1 with the U. S. 
Air Force through the Air Foree Cambridge Research Center, Air Research and Development Command. 
+The principal frequency of a membrane or cutoff frequency of TM modes in a waveguide is de- 
termined by  . 
**The analysis of this section holds only for a convex boundary but may be extended to star-shaped 
domains and perhaps to arbitrary simply connected regions. 
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Translating the origin to any point (h, k) within or on the boundary and introducing the 
new coordinates x’ and y’ gives for the new parametric equations of the boundary: 
x’'(6) = x(@) — h, y'(6) = y(6) — k; 0 < 6 < 2x. (4) 


In terms of new polar coordinates (p, ¢) with (h, k) as pole the boundary is, say, 


p= P(), (5) 
where 
x’(¢) = Pd) cos ¢, y'(¢) = P(¢) sin 9; 0 < @ < 2x. (6) 
Thus, from (4) and (6), 
_, y(0) —k = 
= tan™'' ; ( 
? x(8) — h ‘) 


These simple transformations are shown in Fig. 1. 











We now define the ‘‘a-mean radius with respect to (h, k)”’ of the boundary curve (1) as 
‘ b I . 


ty _ 

r.AR(O))o = {1 | P*(o) do? (8) 
aT Jo J 

From equations (3)—(7) the above definition becomes 


AR). = 4— | (R°(6) — 2R(0)(h cos 6 + ksin 0) + (h? + KY]? 
FT Jo 


. ; ; dR ‘ |“« 
[Ro — R(6)(h cos 6+ ksin @) — 19 (h sin 6 — k cos 0 | do? (9) 
) 
In retaining the limits of integration we have used the periodicity of all functions involved. 
This expression shows explicitly the dependence of ,,.(R(@)). on h, k and R(@). The 
“y-mean radius’, R, , for the boundary curve (1) is defined as follows: 


a > 2, R. = Ming.) v.(R(O))« 
‘gai all (h, k) withinr = R(@). (10) 


a<2 R. = Maxa.s 2.460 8)) «) 
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The points (h, k) which furnish these extreme values will be called the “‘a-centers”’ of the 
cross-section. The case a = 2 is exceptional since ,,,(R(@)). is independent of h and k; 
in fact R, is the ordinary mean square radius or the radius of the circle with area equal to 
that of the given cross-section and any point is a 2-center. For a particular choice of 
(h, k) it can be shown (see [2], paragraph 6.10) that ,,,(R(@)), is a monotonic increasing 
function of a@ [with slight restrictions (see [2], paragraph 6.6 and 6.10) on the function 
R(@), which are assumed to hold]. By means of the definition (10) it is easily proved that 
the same must be true of the a-mean radius, or explicitly, 


Rs he for a<B. (11) 


The equality holds only if the boundary is a circle. In the appendix it is shown that 
R., = afor a circle of radius a and the center of the circle is an a-center. 

Interesting consequences of the definition (10) are obtained if we let a ~+ ©. Thus 
(see [2], paragraph 6.10) from (8) 


,«h(0)). = max, P(¢) | 
for all@ inO < @ < 2z. (12) 
(R(@))_. = min, Pe) 


The corresponding mean radii become, by (10), 


R.. = min max, P(¢)] 
for all (h, k) withinr = R(6) and@inO0 <@ < 2x. ~~ (13) 

iio = Max, min P)) 

The first expression is the radius of the minimum circle containing the given cross-section 

(i.e., the circumscribed circle) and the second is the radius of the maximum circle con- 

tained within the given cross-section (i.e., the inscribed circle). 

3. Bounds for the principal eigenvalue. For a given cross-section (1) consider the 
family of circles with radii r = R, , the a-mean radii of the boundary. The principal 
eigenvalue, \,. , for any circle is the smallest positive number (see [6], section C) such that 
JolXelb« 0, and thus 

ha = ju /Re - (14) 


Here jo; is the first root of the zero-order Bessel function. From the monotonicity of R, , 
expressed in (11), it is obvious that \, must be a monotonic decreasing function of a; 
that is 

: ae for a <8. (15) 


It is a well known theorem that for any cross-section the principal eigenvalue is bounded 
above and below by, respectively, the principal eigenvalues for the inscribed and circum- 
scribed circles. Thus, from (13), this theorem is expressed in our notation as 


he SAS Aves (16) 


Here \ is taken to be the principal eigenvalue of the arbitrary cross-section (1) and the 
equalities hold only if this figure is a circle. Further, Faber [4] and Krahn [5] have proved 
the following theorem: among all cross-sections with the same area, the circle has the 
least principal eigenvalue. By the previous discussion of the case a = 2 this theorem 
becomes in our notation 


Ae SA; (17) 
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the equality again holding only for a circle. From the monotonic property (15) and the 
theorems expressed in (16) and (17) we have for any cross-section 
As SO, for 2 <a < @, (18) 

This may be considered a trivial extension of the theorem of Faber and Krahn to the 
following 

Theorem: Among all cross-sections with the same a-mean radius, foranyain2<a<o, 
the circle has the least principal eigenvalue. 

From the monotonic property (15) it is obvious that for a given boundary the closest 
approximation \, , to the true eigenvalue i, for the above range of a, is dz . If the boun- 
dary is a circle all \,, are equal to the true value, and if not then A, < A. Thus we are led, 
naturally, to seek some a. < 2 such that for any boundary \,. < A. This bound would be 
closer than that obtained from the equal area circle. If such a value could be obtained the 
above theorem could be extended to the range a. < a <@, where a. < 2, in which case 
it would be a non-trivial extension of the theorem of Faber and Krahn. We shall proceed 
to show that for boundaries which are perturbations from a circle such a value does 
exist and we shall obtain it. 

4. Almost-circular boundaries. 

r= R(0) =R, + >d ERO); (19) 


1 


Let the boundary of a domain be given by 


where 

R,(0) = >> |R,; cosné + S,; sin n6], j=1,2,-°-. (20) 
Here « is a small positive parameter which is a measure of the deviation of the curve from 
a circle, and S,; = 0. The problem formulated in (2) has been solved for the above 
boundary, up to second order in e, by various authors (see, [1], Chap. IX, 209-211; 


[3], Paragraph 6.5(a); [6], section C). The principal eigenvalue is found to be 


: (Pp 
Jo1 (Se 

=* l — <— 
. Ry | “\Ro J 


= Ras) ae E J jox) | P _ ' | 
4 (ator) _ hi: lee , Soa 2 4 St, ree 2 
ao (5 R, | 2R2 dX | J J.(jn) + 21%" + ; de (21) 


0 n=1 


Here jo; is defined below (14) and J/(jo.) = (d/dx)J,(x))..;,, . In order to compute A, , 


defined by (14), the a-mean radius of (19) must first be obtained. 

Since the center of a circle is an a-center (see appendix) and the boundary (19) is 
almost circular we assume that its a-center is located at some point, (eh, ek), near the 
“center”. Then using (19) in (9) and expanding in powers of « gives 


on, ex e( 8) a= Ry + €{o,0(R1(8))s} 


} 


+ €)0,0F2(8)): + = 


, ‘ > 2 
QR, [o.o(F21(8)), i oo, (8));] 








et »2" 
+ = fh? +] + ax | teu (—hsin 6 + k cos 0) do 


a2 


— (a — 1)’ | R,(0)(h cos 6 + ksin 6) ao} +¢--- 


70 
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In the first integral above integration by parts gives 


| < R,(0)(—hsin 6+ k cos 0) dé = [ R,(0)(h cos 0 + ksin 6) dé. 


“0 
Now the integrals appearing in (22) can be combined to yield the term: 


(c 2) 
oR, i 


a2r 


R,(6)(h cos 6 + ksin 6) dé. (23) 
Let the first two Fourier coefficients of R,(@) be denoted by 


RR, = i | R,(0) cos 6 dé, 
T Jo 
(24) 


8, = . | R,(6) sin @ de. 


“0 


Using expressions (23) and (24) all terms, up to second order in e, which involve h and k 
in (22) become 
(a — 2) 


iP {h? — 2hR,, + k? — 2kS,,}. (25) 
tllo 


i 


To find R, we must obtain h and k such that the above expression is a minimum for 
a > 2 and a maximum for a < 2. These values are easily found to be: 


h = Ril 


k was 


The a-mean radius of the boundary (19) is thus given by: 


for a>2 and a <= &, (26) 


( 
a R - €} ) R,(0) se oo €'0.06Fa(9)) 








(27) 
a — l) » , 2 ~~ 2 : 
. 2 oR, lo o(Ri(0)), — o.o(?,(8));] - on + Si } + é 
for all a, where the R,, and S,, are given by (24). 
The principal eigenvalue of a circle with the above radius is, from (27) in (14): 
jos o.0(R,(8))s | {[e o(Ri(8)), De | _ o.o(F2(8))s 
Ae = PR, [ es 2° Ro Ro 
(28) 


(a — l) —" , = ; 
_ . [o.0(21(8))1 — 0.0(,(8)) | + oS ”) Re, + sip or oe |. 
Further if each R,(@) is given by a Fourier series as in equation (20) the above principal 
eigenvalue becomes: 


_ dor (Ror\ (Ba)? a 
n=l - de | pt ey(pe Ro — ape Rin + Si) 


(a — |) 2 3 
bat ris > (Ri + sty} + € 


n=2 


(29) 
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To compare this with the true principal eigenvalue for a boundary given by equations 
(19) and (20) we form the difference (A — \..) where J is given by (21); this yields 


2 Jo . . J2( 40 ) a\ >2 v2 \ ¢ 
iat egal EE be ‘jn J Cin.) + 5(@a + Sh) 1+ € --- (30) 
oo 0 m=2 \ ¥ n\ Joi) a) 
Here we have used the fact that 
. J i( Jor) 
sl > (7 ali G1) 


this accounts for the cancellation of the terms involving R,, and S,, . In order that \, be 
a lower bound on the true eigenvalue the difference given by (30) must be positive. Also 
the best such lower bound must be that for which a = a, , the least a for which this 
difference is still non-negative (since \, is a monotonic decreasing function of a). 
To make (30) positive for all R,, and S,, all coefficients of (R;, + Si.) in (30) must be 
non-negative; or 
TnI 1) > a 


Jo1 seaah for ne 2. oy *** fa 
‘ Iii. > 2 ») 


r 


er 6 fr 
a = —2 g.l.b. E sere P 33) 
n=2.3.4,°¢8 J »\ Jor) 


Then ap is given by: 


To find the greatest lower bound in the above expression we use* 


J At 2x" 
‘Fe a- 2s, (34) 


J (x) ’ S Ine coe 


where j,, is the v-th zero of the n-th order Bessel function. Thus we must consider 


‘ Bhd Jor) . 2701 
Wace. ea ZL. a ar n= 2,3,4,--- (35) 
: J (Jor) v=1 Jav Joi 


However, (see [7], p. 508 (2)) for fixed v the j,, are monotonic increasing functions of 7 
’ L4 Jy J g v 


and j,, > jo: forn > 2,v > 1. Thus as » increases each term of the sum in (35), and hence 
the sum, decreases. Therefore the least value of (35) for nm > 2 occurs for n = 2. Using 
this value (33) gives 

J (Jor) 


Oy = —Qjo, A = (4 — 72.) = -1.779 --- 36 
1" Jol Jos) sid " — 


We have thus proved the following result: For any curve in the neighborhood of a 
circle the principal eigenvalue is greater than or equal to that of a circle with the same 
a) mean radius, up to second order in e. This lower bound is the best possible in the sense 
that for any a < a there exist curves for which the inequality does not hold. Furthermore 
this result yields a lower bound which is closer to the true eigenvalue than that given 
by a circle of equal area. These results are expressed in our notation as follows: 


My € As, SA Ae = Se, 
*See [7] p. 498 (1). The argument of this paragraph. is essentially that of [3], paragraph 6.5(b). For a 


different proof see [6!, section F. 
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If an approximation to the principal eigenvalue is desired, rather than a lower bound, a 
value of a Jess than a) should be employed since then the terms in the “error” given by 


(30) will not all be positive, and may cancel. The value a = —2 seems reasonable, for 
too small a value would yield a negative error. Furthermore, the calculation of the 
a-mean radius is relatively simple for a = —2. 


We further conclude that if there exists an a, < 2 which furnishes a lower bound for 

all cross-sections, not only those near a circle, it must lie in the range 

ay = (4 — joi) S a < 2. 
If not, the above example shows that we could assign Fourier coefficients (and hence 
determine a boundary) in (30) such that A... > A. 

5. Application to rectangles. The eigenvalues for 1 rectangular cross-section are 
known exactly and thus provide a test of the approximations just obtained. 

Let the rectangle be located with its axes of symmetry coincident with the coordinate 
axis and its sides of length 2 and 20 respectively. Then the boundary may be expressed 
in polar form as r = R(@) with: 

] 


2/ —_ ete f 
R(6) cos 6” B< 6<P, 
if oC 

Hp) = ——, B<6<7-8B 

sin 6 

(37) 

R(6) = oat r—-B<6<74+8, 

cos 6 
R(@)=-—~, wrt+B<0< 2-8, 

sin 6 


where 8 = tan™ go. 
The a-mean radius becomes, from (9), (10) and (37): 


R, = 2 | Le+ef ao). (38) 


0 cos* 6 /s  gin* 6 





Here we have used the assumption that the a-center lies on a line of symmetry if such 


exists. The a-mean radii are, for a = 2, 1, —1, and —2: 


1/2 
n= (2) 
T 


— 21/2 _ (1 + 0°)" 74] 
r, =2|ma +e) eli ya) 
(39) 
> 7 TT o 7 
ae = 2(1 Ys ee 


R.. = ()""4| 6 + Bo? + (z - yl 


From these we may compute the principal eigenvalues of the circles with the corre- 
sponding radii (by means of (14)). The exact principal eigenvalue for the rectangle is 
(see [3], P. 253) 

x (1 +- o 1/2 

=_ |e ee atl (40) 


2 o 


A = 
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We thus obtain the following ratios: 


Xo l o 
r 1 ( "T ' oe ayiss — 
= 1 + In (1 z = —— 355 ——_ 
X, 7m n(l+oa) olny ney |) 
(41) 
ole 
MN. 34 


r ” la 1/2 : eo - 
A = F— (1 + 0°)""[o + Bo" + (x/2 — BI; 8 = tan™ o. 
A2 “Jo 


These ratios give a measure of the closeness of the approximation by various circles as a 
function of the ratio of the sides of the rectangle (see Fig. 2). It is interesting to note that 
the (—1)-mean radius has a constant ratio of approximation for all rectangles. This is 
all the more remarkable since Polya and Szegé (see [3], paragraph 5.9(b)) have obtained, 
for a rectangle, an upper bound on the principal frequency (eigenvalue) which also has a 


constant ratio of approximation. This ratio is 


\* 7 ( aye) — x 
» Ny it 6c 


where \* is the upper bound on the principal frequency. This value is exactly that of the 


asymptote to, the \/A_2 curve of Fig. 2. 





x 
Aah 
1.50; 
a=2 
1.40/- (Equal area appx.) 





0.90-- 





f'~ 

we 
b 
a 
q 





Fic. 2. Ratio of Approximation 
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APPENDIX 


Mean radius of acircle. In finding the a-mean radius of an almost circular boundary 
it has been assumed that the a-center is in the neighborhood of the center of the unper- 
turbed circle. This assumption is valid if it can be shown that for a circular boundary the 
center is an a-center. It is easily shown that the center furnishes a relative maximum 
(a < 2) ora minimum (a > 2) but from the definition given by equation (10) we require 
an absolute extreme. It is shown below that the center of a circle is an a-center for all a. 

By symmetry we may place the center on the z-axis when the boundary is a circle 


given by R(@) = a (see Fig. 3). This yields for the “‘a-mean radius with respect to h” 


yA 











Fic. 3 


+ [ Pw) aa} ) 
(A.1) 
- { [ [a? — 2ah cos 6 + h?|*~*”|a? — ah cos 6) de? , 


\Qr /0 


where —a < h < a. If we let a—+~, which may be done for all h since P,(¢) vanishes 
at one point at most (see [2], paragraph 6.6 and 6.10), there results: 


a+ lhl, 


,(a). = maximum P,(¢) 
(A.2) 
g=— | h a 


(2) 


minimum P,(¢) 
eo 


Furthermore we have by the remarks under equations (10): 
n(Q)o = da. (A.3) 


Using equations (A.2) and (A.3) and the fact that ,(a), is a monotonic increasing function 


of a we must have: 


asa). at {hi, a>2 
—a<h<a. (A.4) 


a—|h|< a). <a, a<2 
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The a-mean radius defined by (10) is now from (A.4), 


R,. = min, ,(@). = 4; a. >, 
(A.5) 


2. = max, ,(@). = @; e< 3. 


This value is attained in either case for h = 0; i.e., the center of the circle. 


REFERENCES 


. Rayleigh, Theory of sound, Vol. I, Dover, 1945. 

. G. H. Hardy, J. E. Littlewood and G. Polya, Inequalities, Camb. Univ. Press, 1934. 

. G. Polya, and G. Szegé, Isoperimetric inequalities in mathematical physics, Princeton U. Press, 1951. 

. G. Faber, Beweis dass unter allen homogenen Membranen von gleicher Flache und gleicher Spannung die 
kreisformige den tiefsten Grundton gibt, Sitzungsberiche, Bayr. Ak. Wiss, 1923, 169-172. 

. E. Krahn, Uber eine von Rayleigh formulierte Minimaleigenschaft des Kreises, Math. 94, 97-100 (1924). 

. H. B. Keller and J. B. Keller, Eigenvalues of nearly circular wave-guides, Res. Rpt. TW-17, NYU, 
Math. Res. Grp., Wash. Sq. College of Arts and Sci., 1952. 

7. G. N. Watson, A treatise on the theory of Bessel functions (Second edition), Camb. Univ. Press, 1948. 


mm whD 


a o 








151 


. 1 FECTS OF SURFACE TENSION AND VISCOSITY ON TAYLOR INSTABILITY* 


BY 
RICHARD BELLMAN AND RALPH H. PENNINGTON 
Rand Corporation, Santa Monica, Cal. 


Abstract. ‘The model used is that of two fluids of infinite depth, with the interface 
initially in the form of a sine wave with amplitude small compared to wave length. 
The fluids are considered incompressible, and only the linear terms in the equations of 
hydrodynamics are used. The first four sections discuss the effects of surface tension 
and viscosity. The fifth gives a few numerical results to illustrate the main points of 





the preceding sections. 

Introduction. If two different fluids having a common plane boundary are accel- 
erated in a direction perpendicular to the boundary, any small irregularities in the 
boundary will tend to change in shape. If the acceleration is directed from the more 
dense to the less dense medium, the irregularities will tend to smooth out (in the absence 
of external forces). Thus the plane configuration of the interface is a stable one. This 
can be illustrated by the usual example of a glass of water sitting at rest. If one considers 
the force of gravity to be replaced by an acceleration which produces the same effect, 
the water and the air are undergoing an upward acceleration. Since the acceleration is 
from the more dense to the less dense medium, the air-water interface is stable. 

Returning to the general case, if the acceleration is directed from the less dense to 
the more dense medium, irregularities of the interface will tend to grow. This is the 
effect known as Taylor instability. An example is the case of glass of water turned 
upside down. Here again the force of gravity may be considered to be replaced by an 
upward acceleration. The acceleration is from the air to the water, and the air-water 
interface is instable. The water, instead of maintaining a nearly plane lower surface as 
it falls, will tend to jet out into long spikes. It is the formation and rate of growth of 
these spikes which we wish to investigate, taking into account the effects of surface 





tension and viscosity. 

1. Taylor’s results. Let us begin by presenting an account of the work done by 
Taylor himself [4] (see also [2]). The model used is that of two fluids of infinite depth. 
The interface (neglecting perturbations) is the plane y = 0, the y axis being vertical. 
The initial perturbation will be of the form cos kx, with amplitude small compared to 
wave length. The problem is then two-dimensional, and the true equation of the interface 
at any time is y = n(2, t), where the function n(z, ¢) is to be determined from hydro- 
dynamic considerations. The fluids will be considered to be incompressible, and only the 
linear terms in the equations of hydrodynamics will be used. 

The linearized hydrodynamical equations in either fluid are 


u, +v, = 0, (1.1) 

U, + . p. = 0, (1.2) 
p 

(1.3) 


1 
utontgtn = 0. 
*Received June 8, 
residence at Princeton University. 


1953. The results in this paper were obtained in 1951 while both authors were in 
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Here, as usual, uw and v denote the components of velocity in the x and y directions 
respectively, p the pressure, p the density, g the acceleration of gravity, and g, the 
upward acceleration of the system. These equations have solutions of the form 


u= —®¢., v= —¢,, (1.4) 
P=P—-(9+ nr + , (1.5) 
where ¢,, + ¢,, = 0 and pp is the mean pressure at the interface in the unperturbed 
condition. 
For the upper fluid we take 
¢@, = Ae *"f(t) cos Kz, (1.6) 
Pi = Do — (9 + 9i)Piy + pili): (1.7) 
and for the lower fluid, 
¢. = —Ae*"f(t) cos Kz, (1.8) 
™~ = (g + 91) Poy + P2(do) ’ (1.9) 
the lower fluid being the more dense, i.e., p2 > p, . The above relations satisfy the con- 
ditions that velocities are finite at y = © andy = —o, and that v, = v, at the (approxi- 
mate) interface. 
The free boundary condition at y = n(2, ft) is 
v—nu—7, = 0 (1.10) 
which upon neglecting the non-linear term yields n, = v = KA/f(t) cos Kz, or 
at 
7 = KAl | f(t) at) cos Kx. (1.11) 
J ts 
The pressures at the interface must satisfy the condition p, = p, . Substituting 


from (1.7) and (1.9) we obtain, after some slight manipulation, 
—(g + g:)(p2 — pi) K f(t) — (a2 + a) f’() = 9, (1.12) 


so that we may take f(t) = sinh nt, since this choice makes the fluid velocity zero at 
t = 0. The value for n° is given by 


' —( ipo — pK 
e (g + gi)(p2 — pi) 1.13) 
Po + Pi 
and the interface is given by 
n = KAn™ cosh nt cos Kx (1.14) 


If (g + g,) is negative, there is a positive value for n. The disturbance grows like 
cosh nt, which means that the motion of the interface is instable. This instability exists 
for all positive K, i.e., for all wave lengths of the initial disturbance. Note that the smaller 
the wave length (= 27/K), the more rapid the growth of the disturbance. This limits 
the use of the above result for arbitrary disturbances. Since the differential equations 
used were linear, one would expect to discuss an arbitrary disturbance by Fourier 
analysis. Let the surface at time t = 0 have the equation 


y = f(z) = > ax cos Kz. 


K=0 
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Then at time t, we have 


y = >> ax cosh nt cos Kz, 
K=0 
where n is essentially K’’?. For t # 0, the series will diverge unless the convergence of 
>" %-0 ax cos Kx is extremely rapid, since cos [h(Kt)'””] grows so rapidly. 

2. Viscosity. The effects of viscosity on the arguments of Section 1 are clear intui- 
tively. Viscosity is not to be expected to remove the instability, but only to reduce the 
rate of growth of the amplitude of the disturbance for any particular frequency. The 
amount of this reduction for small wave lengths is rather startling, however. In particular, 
as the frequency — @, the rate of growth of amplitude — 0. 

The model to be used here is that of Section 1. The (linearized) equations governing 
the motion of an incompressible, viscous fluid are 


u.+v, = 0, (2.1) 
u, + ew =- (uzz + Uy), (2.2) 
p p 
1 m 9 
Ve ~~ Py + g + Ban (v.. + Uyy)» (2.3) 
p p 
where u is the coefficient of viscosity. 
These equations are satisfied by 
u= —¢, — Wy ’ , = —d, + Vz ’ (2.4) 
p= Dp — (g+ ge, + ed): , (2.5) 
provided that (ef. [1, 2]) 
b+ bm =0, 7 (Yer + Ww) = We (2.6) 
For the upper fluid we take 
¢, = Ae~*"*™ cos Kz, (2.7) 
y, = Be™**™ sin Kz, (2.8) 
Pi = Po — (9 + gi) + rid): (2.9) 


where 


9 *? nn 
at = Ke +f. 
Mi 


For the lower fluid, 


o. = Ce*"*™ cos Kx (2.10) 
v2 = De™”*™ sin Kz 2.11) 
P2 = Po — (9 + Gi)Pay + polde): (2.12) 


where 


2 r2 P2n 
Mm = + 
Me 
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In order that the velocity components remain finite for y — +, it is necessary that 
the real parts of m, and m, be positive. 
Let the interface be given by y = n(z, t). Then n, = v, or n, = K(A + B)e™ cos Kz, 


from which we obtain 


n = K(A + B)n™'e™ cos Kx (2.13) 
The boundary conditions at the interface are 
us = Us, Vv; = V2 (2.14) 
Ov, Ov. = 
—D, + 2p; ay = —P2o + 2po ay (2.15) 
Ov Ou, Wo OU2 
(2 + a) ~ po( 2 ” dus) (2.16) 


The last two equations state the equality of the components of the stress-tensor. Sub- 
stitution in Eqs. (2.14) to (2.16) gives four conditions on the constants A, B, C, and D. 
From (2.14) we obtain 
KA +m,B —- DC+mD=0 (2.17) 
A+B+C-D=0 


From (2.1: 


anil ab "] _— VK = —(¢ g1)( Nahas K , 
| go Ps e = oie ou K? | A +" |= 9 F g\os = ps) - — 2u,K m, [B 


)) we have, after some simplification, 


fi i. if) 
+ [pon + 2u.K°|C — 2u.Km.D = 0 (2.18) 
Henceforth, we will let —(g + g:)(p2 — p,)K = 8. Similarly, using (2.16) and some 
algebra, we obtain 
Qu,K°A + w(K? + m)B + 2u.K°C — p.(K* + m)D = 0 (2.19) 
These four equations are linear and homogeneous in A, B, C, and D. They have non- 
trivial solutions if and only if the determinant of the coefficient vanishes, 
l l | 1] 
| ~ = | 
| K mM, —K Me 
| | 
| 1=0 (2.20) 
| Pu, K- u,(K? + mi) Qu.K" —p(K* + m3) 
B 9 1 B 9 K 2» K? 9 K 
—~ — pin — “pith — 26, AM, Pot - “fe — Sh2ik Me 
Ln” n = 


This equation reduces to 
[—8 + (pi + p2)n |[(uiK + [2 M2) a (uk > Him;)] 
+ 4nK[uwiK + pwome)[u.K + w,m,] = 0 (2.21) 


Equation (2.21), using the values of m, and m, given above, yields an equation of tenth 
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degree in n. Since the roots cannot be directly determined, it will be more profitable to 
avoid rationalization and see what information can be obtained by other means. 

In Section 1 we found that n was positive, which implied that instability occurred, 
when (g + g,) was negative. The value of n which determined chiefly how fast the 
amplitude of the disturbance grew was 


7\ 1/2 
oe 4(=2+ ar ~ eK) (2.22) 
Pi + po 
In the present case, then we expect that when 8 = —(g + g,)(p2 — p,)K is positive, 


there will be at least one root of equation (2.21) with positive real part, and that for 
this root, 


Re(n) < ( = | .) : -_ (=o + 2X — dK)" (2.23) 
(Pp: + po) P2 + py 

We shall determine whether equation (2.21) has a root with positive real part by 
considering n a complex variable and applying the principle of the argument to the right 
half plane. So long as n remains in the right half plane, the quantities m, = 
(K*® + p.n/u,)'? and m, = (K* + p.n/p,)'” remain on one branch of their domain of 
values, so we will run into no confusion if we write equation (2.21) as 


[—8 + (p1 + p2)n [mK + (usK* + open)” + mK + (ui K* + ipin)'”*] 
+ 4nK[uiK + (u2K? + popon)'*|[u2K + (uiK? + uipin)'”*] = 0 (2.24) 


Let us use the contour C consisting of the part of the imaginary axis between (0, R) 
and (0, —R) and the semi-circle in the right half plane with this segment as diameter. 
If we denote the left-hand side of equation (2.24) by f(n), the principle of the argu- 


ment states that 


* f’( , E : . 
| 7 n dn = i[change in argument of f(n) around C] 
Jo f(n 
(2.25) 
= 2ri{number of zeros of f(n) within C] 
so that 
; , f'(n) : ——" 
lim | L dn = 2ri{number of zeros of f(n) in right-half plane] (2.26) 
rao vc J(n) 
Now 
° fl ( ) ; ax/2 , ,1¢ . : , -iR , 1 “ 
lim | LY dn = lim | f'(Re") i Re’* de + lim [ f'n) dn (2.27) 
paw do JW Row /— 9/2 f(Re’*) raodirn Jn) 


provided that the limits on the right-hand side exist. Since they do, we shall evaluate 
them separately. To evaluate the first term, note that the highest power of n appearing 
is n°”. In the limit this is the only term which will matter, so 


»*/2 ¢r/ PD, i¢ n/2 Jy, 38/2 


, sw ee : ae se 
lim | — i Re’® dé = lim —’ Re** do 
Roo /-2/2 (Re'®) Row /-—#/2 n” 


/ 


ee/2 = Ip ,16)\ 3/2 ep w/2 
= lim [ 2A) Rede = [| 5/2idd = Swi/2. (2.28) 
Row J— 4/2 (Re**)*” J—#/2 
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Hence, in the limit the change of argument of f(n) when n traverses the semicircle 
is 54/2. The second term in equation (2.27) is 7 times the change of argument of f(n) 
as n traverses the imaginary axis from +io to —io. This change of argument can 
be seen directly. 

We next ask for an upper limit on the value of the real part of this root. We note 
first of all that the root itself is real. For positive real n, f(n) is real and continuous in 


mn for all K. For n = 0, f(n) = —28(u, + u.)K <0. Forn = (8/(p: + p2))'””, f(n) = 
4nK[uiK + (u3K? + popon)*|[u2K + (uiK° + uipin)'””] > 0. Hence there is a positive 
real root between n = 0 and n = (8/(p; + p.))'””. We already have, then, an upper 


limit on the root. At the root, n < (8/(p; + p2))'””. This is just as expected, for the 
positive value of n when viscosity was neglected, see (1.13), was (8/(p: + p2))'”*. 

We can find an upper bound on this root which shows more about the nature of the 
root. To do this, we rewrite (2.24) as 


1 
\ ( — \K 4 
(gi: + 9)(e: pi)K + (pr + pa)n” ‘| iK +} (kK? + b2psn)' 


a —_—= ' so 999 
3 Mok + (u2?K? + pp,n)” | Tie = aa 


keeping in mind that g + g, is negative. Consider also the comparison equation 


(9. + g)(p2 — pi)K + (po, + eae'l| a a | +42K=0 (2.30) 
mK + wK mK + mK 
For any K, the positive root z of (2.30) must be greater than the positive root n 
of (2.29). The second factor of the first term has been increased, and at a root this must 
be counterbalanced by an increase in the second term or a decrease in the first factor 
of the first term, both of which require an increase in the root. For an upper bound 
on n, then, we have only to give the value of z. We rewrite (2.30) as 


(pr: + poe” + 22K7(ur + we) + (9 + gi(er — aK = 0 (2.31) 

The positive root is 
lus + oe) K* + (Cus + 2)*K* — (9 + gi)(o2 — pi)(02 + 01)K)™” (2.32) 
Po + Py Fae 


The most interesting thing about this root is the fact that it has a maximum for 
some K. Thus the introduction of viscosity has eliminated the tendency for disturbances 
of small wave-length to increase without bound. We would like to know the value of z 
at this maximum, since this value will be an upper bound on n for all K. Differentiating 
(2.32) with respect to K and setting dz/dK = 0, we obtain 


42(u, + be) 
Substituting this value in (2.32), the result i 
[-(g ats (p2 — pid)” 
(2.34 
2p: + 178 (uy + ue)” ) 
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This occurs for 


- . 1/3 ie 
- ola anle = ell fot po (2.35) 





K 
so that for all K we have 


n < Cet ale — ool (2.36) 
2(p1 + po) (ur + M2) 

One can, of course, make better approximations for n.- For the general case this 
process does not seem to offer much, since the general state of affairs is now established. 

A quite complicated, but straightforward, calculation shows that n has only one 
maximum. 

Note that in this present Section one cannot satisfy the condition that the velocities 
be zero when ¢ = 0. Apparently because of the linearization performed, one obtains 
no motion at all if one attempts to satisfy this condition. 

3. Surface tension. We now introduce the effects of surface tension into the argu- 
ments of Section 1, ef. [2]. It is to be expected that the presence of surface tension will 
remove the instability for sufficiently small wave lengths. This is indeed the case, as 
will be shown. 

To introduce surface tension into the arguments of Section 1, we have merely to 


replace the condition p, = p, by 
P2 - Pr + Tin = 0 (3.1) 
Substituting as in Section 1, we obtain 
—(g + gi)(p2 — piln + polde)e — pibr): + Tinee = 0 (3.2) 
or 
—(q + 9:)(p2 — p;)AKn™ sinh nt cos Kx — (p2 + p,)Ansinh nt cos Kx 
— T,AK*n™' sinh nt cos Kx = 0, (3.3) 
so that 
»_ —(g + gid(p2 — pi) p T , 
ee 3.4 


The condition given for Taylor instability was that g + g, be negative. But we see 
from equation (3.4) that the amplitude of the initial disturbance grows only when 








—(9 + goo — or)  _ _TiK* r 
P2 + pi * Pi + pro iis G5 
or 
1/2 
K < [e+ ale — 09) 3.6) 
or 
> 2r/ P, 7 (3.7) 
m —(g + 9:)(P2 sia pi) , 
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where \ = 27/K is the wave-length of the initial disturbance. Thus for wave-lengths 
smaller than those satisfying condition (3.7), there is no instability.* 

Another fact of importance is expressed by equation (3.4). Since the right-hand 
side has an absolute maximum, there is a “most dangerous frequency,”’ i.e., a frequency 
for which the amplitude of the disturbance grows most rapidly. 

The most dangerous frequency is that frequency for which n, or n’, is a maximum. 


At this frequency, then, 


d (g + 91) (po — ~) = T, P| 90 
oa — a — = oo (; Ss 
dK | Po + pi . Po + pi ’ si 


from which 


1/2 


[—(g: + g)(e2 — pi)] 


K= TE (3.9 
a ) 
Substituting this value in equation (4), we have 
2 a 9:)(p2 — p,)]’* 
ee [—(g + g:)(p2 — pi)] (3.10) 


3(37,)'” Po + pi 


It is remarkable to note the small effect which the numerical value of the surface 
tension has on the rate of growth of amplitude. Although it is the quantity which places 
a limit on the rate of growth of amplitude, it is felt numerically only in the one-fourth 
power, as equation (3.10) shows. 

4. Viscosity and surface tension. In this Section we combine the results of the two 
preceding to give an over-all picture including both surface tension and viscosity. We 
would expect that as in Section 3, there would be no instability for small wave lengths; 
and that for longer wave lengths, the rate of growth of amplitude of the disturbance 
will be less than that given in Section 3. 

The procedure will be to take the arguments of Section 2, where viscosity is considered, 
and alter them to include the effects of surface tension. To do this we must replace 
equation (2.15) of Section 2 by 


OW» Ov i 3 
=p, + Qe + p: — 2 —* — T, 4 = 0 (4.1) 
oy Oy Ox 
Substitution in this equation as in Section 2 yields 
a, o a 4 
e —+ pn + 2u,K a 4. E : +t 2uKm, |B 
n D 
+ [—2u.K° — pn|JC + 2u,.Km,D = 0, (4.2 


where a = p(g + g:)(p2 — p,)K — TK’. 


*This explains the hanging of water droplets on the underside of a horizontal surface, such as a ceiling. 
Such a droplet is undergoing an upward acceleration of 980 cm/sec? and will tend to drip because of 
Taylor instability unless its effective wave-length is too small to satisfy (3.7). For water, the critical 
wave-length is about \ = 2x(74/890)!/2 = 1.73 cm. Droplets of larger diameter will tend to drip, while 
smaller ones will tend to hang. (Actually, of course, the true critical diameter will be different because of 
circular symmetry, etc., but the above at least contains the principle involved.) 
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The other three conditions on A, B, C, and D are the same as those in Section 2, namely: 
A+,B+C-D=0 
KA + ™m,B — KC + mD = 0 (4.3) 
Qu,K°A + w(K? + m)B + 2u.K°C — w(K? + m)D =0 


Equations (4.3) are linear and homogeneous in A, B, C, and D. They have non- 
trivial solutions if and only if the determinant of the coefficients vanishes, 


! 1 1 1 | 
| 
K m, —K Me 
| =0 (4.4) 
Qu,K? u,(K* + mi) Quek? —p,(K? + m3) | 
| 
—~ +n +%K* — : + 2u,Km, —2yu.K* — pon 2u.K mp | 
il od 


This equation reduces to 
[—a + (p, + po)n™)[(uiK + wome) + (uoK + wm)] 
+ 4nK[wiK + pome)[u2kK + wim,] = 0 (4.5) 


This is precisely equation (2.21) of Section 2 with 8 = —(g + g;)(p2 — p,)K replaced 
by a = —(g9 + g:)(p2 — pi)K — T,K”’. In Section 3, where surface tension alone was 
considered, we found stability for a < 0 and instability for a > 0. We shall show that 
these conditions still hold. 

For a > 0, the result is immediate from Section 2. Paraphrasing the results of Sec- 
tion 2 for a > 0 instead of 8 > 0, we may state: for a > 0, equation (4.5), where n,m = 
(uiK° + wipin)'* and pom, = (u3K° + popon)'””, has just one root with positive real 
part. This root is real and less than (a/(p; + p,)'””. We will return to the problem of a 
better estimate of this root after proving stability for a < 0. 

To establish stability for a < 0, we again apply the principle of the argument, as 
in Section 2. The result is established by a series of straightforward but laborious argu- 
ments which we shall omit. 

We turn now to the instable case, a > 0. We found that in this case the equation 


, ee oe: ee ee One 
ae as ae l= + GiK? + npn)? + aK + Gk? + | 
+ 4nK = 0 (4.6) 


had one positive root. There are two immediate upper bounds for this root. The first, 


already given, is 














a id (= + 9:)( po sss p,)K — ss - 
LS = co ne oe hoe 4. 
nt Pi + peo (4.7) 
and for all K, 
a [ g+ 9:)(p2 tek p,)]°”* 
n< 3i/4 PS, +p)? a (4.8) 
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Relations (4.7) and (4.8) state merely that the rate of growth when both viscosity 
and surface tension are considered is less than that when surface tension alone is con- 
sidered. 

The second upper bound on n comes from comparison of (4.6) with equation (2.30) 
of Section 2. Since 

a = —(9 + g:)(p2 — p)K — T,K* < (9g + g:)(p2 — pK, (4.9) 
the root of (4.6) must be less than that of equation (2.30) for given K. This may be 
seen in the following way. Suppose the value of a@ in (4.6) is increased. The first factor 
of the first term tends to become more negative. An increase in n will decrease both 
factors of the first term and increase the second term, to counter-balance the change in a. 
Thus the root of equation (2.30) is an upper bound on the root of (4.6). This is merely 
a statement of the physical fact that the rate of growth when both viscosity and surface 
tension are considered is less than when viscosity alone is considered. 

From the study made of equation (2.30) we can give an upper bound for the root of 


(4.6), namely 


at + pe) K* + ((ur + ue)*K* — (g + gi)(O2 — pr)(o2 + p:)K)'”” (4.10) 
(p2 + pi) ; 
and for all A, 
—(g 4 (po — we 
n < Lg + g)(p2 — 11)! “1p 


2(p1 + po)” *(ur + be)’”” 

The upper bounds on n given by (4.7)-(4.10) will not usually be of great practical 
value. For particular cases, numerical methods must be used. 

A little can be said about the frequency for which (4.6) has a maximum root. The 
effect of viscosity is to shift the maximum towards smaller K, or greater wave lengths. 
Furthermore, n has a unique maximum as a function of K. 

5. Numerical examples. In order to demonstrate the effects of surface tension and 
viscosity, we give some examples for ordinary fluids. 

Example 1. If the two fluids involved are air and water, surface tension would be 
expected to play an important role in the development of Taylor instability. We use 
0, Pearer = lg/ec, T, = 74 dynes/em, g + g, = —2 X 10* cm/sec’ = —20g. 


water 


Pair 
Figure 1 shows values on n vs. k when surface tension is considered and when it is 
neglected. The corresponding equations are n° = 2-10 K — 74K* and, n° = 2 X 10°K. 

For the surface tension case, n has a maximum of about 355 at K = 9.5 (A = 0.66 cm) 
and drops to zero at K = 16.4(A = 0.38 cm). The deviation from the no-surface tension 
case is indistinguishable for K < 3(A > 2.1 em). 

Experiments have been made by Lewis [3] for accelerations on the order of that 
used above, at wave-lengths on the order of one centimeter. However, the published 
results are not in a form which allow comparison with those given above. It would appear 
that experimental verification of the effects of surface tension should not be difficult to 
obtain with apparatus like that used by Lewis. 

Example 2. If the two fluids involved are air and glycerine, both surface tension 
and viscosity would be expected to play an important role in the development of Taylor 
instability. We use 

Pair = 0, Peiycerine = 1.26 g/cc, Meir = O, Meiycerine = 14.9 poises, 


T, = 63 dynes/cm, g+ 9: = 2 X 10‘ cm/sec’. 
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Figure 2 shows values of n vs. k under four different conditions: 


1. Neither surface tension nor viscosity acting. 
2. Viscosity ol ly acting. 


3. Surface tension only acting. 
4. Both viscosity and surface tension acting. 
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In this way the relative importance of the two effects for various wave lengths are 


made apparent. The corresponding equations are: 


1. n>? =2-10°K 


; - 9 l l 
9 = ( oR) i. 9Ryn* = cio soe ~ 2 
a. [-3 + USHA + 1 I 14.9)"K* + (14.9)(1.26)n)'” hx | 
sul ink = () 


# “ee os is 
3. n =2-10K — 1.96 


. ; ” 
_ [—2+ 101.26) 63K? + 1.26n?]| -——=>— nn ' 
. L.2O)K + 65K ' | aoe + (14.9)(1.26)n)'? * x! 


It is seen that the viscosity is unimportant for K < 1(A > 6.28 em) and that the surface 
tension is unimportant for K < 3(\ > 2.1 em). Experiments have been made by Lewis 
[3] for accelerations on the order of that used above, at wave-lengths on the order of 
one centimeter. It would seem that the viscosity effects would be apparent in these 
experiments. This would lead to an observed value of » much smaller than that pre- 
ducted by the theory for non-viscous fluids. However, the experiments gave an observed 
value n greater than that predicted by the simple theory. Lewis explains this on the 


basis of viscous drag on the channel sides in the apparatus. 
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LATERAL VIBRATIONS OF TWISTED RODS* 
By 
ANDREAS TROESCH (New York University), 
MAX ANLIKER ann HANS ZIEGLER (Eidg. Techn. Hochschule, Ziirich) 


1. Introduction. Figure 1 shows part of a straight rod which is twisted, like a 














Fic. 1. Twisted rod. 


propeller blade, in the unloaded state. The bending vibrations of such a rod have been 
dealt with by H. Reissner,’ who derived the differential equations under most general 
assumptions, without, however, to the best of our knowledge, proceeding to their solution. 
Subsequently, a number of authors have treated particular cases by approximate or 
experimental methods.’ In view of the possibilities offered by modern computers, the 
problem actually is ripe for rigorous treatment. It is solved in this paper for an isotropic, 

*Received Sept. 11, 1953. 
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homogeneous rod built in at one end, under the assumptions that the mass per unit 
length, u, the twist per unit length, w, and the principal flexural rigidities, a, 8, be constant 
(one of them being taken infinitely large). 

The eigenvalue problem, which is easily established by a method recently used in a 
stability investigation,” is of 8th order. The eigenvalue equation is obtained by equating 
an 8-row determinant to zero, the elements of which contain the eigenvalue in the form 
offan exponent. For small and large values of the total twist, ¢, the natural frequencies 
of the rod can be computed by development. In the intermediate range, they have been 
evaluated by means of the sequence controlled computer of the Institute for Applied 
Mathematies of the Eidg. Technische Hochschule.* 

2. Reference frames. Figure 2 shows the deflection curve of an originally vertical 


z Ww § 








x 
Fic. 2. Coordinate systems. 
rod, S denoting the center of gravity of a generic section given by the are s. The curve 
can be referred to three different coordinate frames, viz., 
(i) the fixed frame z, y, z, defined by the principal axes of an end section and the 
axis of the straight rod; 
(ii) the principal frame &, 7, ¢, 
together with the tangent of the deflection curve; 
w with vertical axis w, obtained from &, 7, ¢ by rotation 


defined by the principal axes of the section at s, 


(iii) the raised frame uw, »v, 
about the nodal line, i.e. the horizontal line in the plane £, 7 passing through S. 

Let t denote the unit vector along the tangent of the deflection curve, while a is 
an arbitrary vector the first two components of which are small. Provided that the 
inclination of the deflection curve is small, the components of a in the last two reference 
frames are connected by the relations” 

a, = a; + ta; , a, = a, + t,a; , ay = ay. (2.1) 
3H. Ziegler, ZAMP 2, 268 (1951). 
‘The authors are indebted to the Director of this Institute, Prof. Dr. E. Stiefel, who put the com- 
puter at their disposal, and to Prof. Dr. H. Pallmann, President of the Schweiz. Schulrat, for granting 
ivourable conditions. Also, they are obliged to Dr. H. Rutishauser for valuable advice with 





them very f 
regard to the establishment of the computing program. 
5These follow at once from Fig. 2; see also H. Ziegler, loc. cit., (4.1). 
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If, for the present, s is interpreted as a measure of time, the principal frame slides 
along the deflection curve with unit speed. Its state of motion consists of a translation 
of velocity t and of a rotation of angular velocity x + wt, provided that w = dy/ds 
denotes the angle of twist per unit length of the unloaded rod, while the components 
of the deformation vector «(xz , «, , k;) represent in turn the curvatures of the deflection 
curve [measured in the principal planes n, ¢ and &, ¢] and the elastic twist per unit length. 
In the raised system the components of x are, according to (2.1), 


K, = Ke + tykz, K, = Ky + t,ky , Ky = Kr, (2.2) 
hence, x + wt has the components 

ke + (ky +o), k, + t,(x; +), ky +. 
The first two components describe the motion of the principal frame with respect to 
the raised one, whereas the last component represents the rotation of the raised frame. 

3. Differential equations of motion. Let k denote the external force and m the 

external moment, both taken per unit length. If K and M represent the internal forces 
in section s, the conditions of equilibrium for an element of unit length are 


A +k=0, MitxK+m=0. (3.1)° 


Since the angular velocity of the raised system is u(0, 0, x, + w), the equilibrium con- 
ditions, for an observer taking part in this motion, are, instead of (3.1), 


dK 


 tuxK+k=0, Mm tuxMt+txK+m=0. (3.2) 


Resolving them and denoting derivatives with respect to s by primes, we obtain 


Ki — (x, +o)K, +k, = 0, 


Ki + (x. + w)K, + k, = 0, (3.3) 
KL +k, = 0; 
M! — (x, +o)M, + t,.K, — t.K, + m, = 0, 
M! + (x, +)M, + t.K, — 4.Ky +m, = 0, (3.4) 
M: +iK, — tA, +m, = 0. 


Denoting by a, 8, y the flexural rigidities with respect to the principal axes and the 
torsional rigidity, we have further, provided that the strains remain small, 


M; =ak:, M, = Bx, ; M; = yx; . (3.5) 


When the rod executes free vibrations, k is the inertia force per unit length and m 
its moment. For small vibrations, k, , k, and m, can be treated as small quantities of 
the first order. Excluding longitudinal vibrations and the higher modes of flexural 
vibrations, we are justified in replacing s by z and in suppressing the quantities 
k. , m,, m, , since they are small of higher order. So is K,, , while K, , K,, M., M,, 
M,, , tu, ts, Ku Key K» are of the first order and t, = 1. 


6See, for instance, H. Ziegler, loc. cit., (2.1). 
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If terms of the first order alone are retained, the third equation (3.3) becomes trivial. 
The third equation (3.4) reduces to Mi, + m, = 0, thus, torsional vibrations can be 
treated separately. Excluding them, we have m, = 0, M, = 0. From (2.1) follows 
m,; = 0, M; = M,, M, = M,, M; = 0 and, owing to the third equation (3.5), 
Ke = Ku, Ky = Koy Kr = Kk» = 0. (3.6) 
Thus, the four remaining equations (3.3), (3.4) reduce to 
Ki — oK, + k, = 0, Mi — wM, — K, = 0, 
(3.7) (3.8) 
Ki + oK, + k, = 0, Mi + oM, + K, = 0, 


while the first two relations (3.5) may be replaced by 
M, = ak, , M, = Bx, . (3.9) 
According to Section 2, the angular velocity u* of the principal frame with respect 
to the raised one has the components 
Ke + tu( ke a w), Ky + t,( ky + w), 0. 
Owing to (3.6), it can be written u*(x, + wt, , x, + wt, , 0). For an observer moving with 
the raised frame 


dt 


= u* X t; 
ds 
hence 
.=-—wt, «=t—awt,. (3.10) 
If, finally, r denotes the radius vector of S (Fig. 2), 


j 
t=o +u Xz; 
aus 


thus 
, =" — o;., t, =ri+or,. (3.11) 


u 


Eliminating x, , x, and t, , t, from (3.9) to (3.11), we obtain the bending moments 


M, = —alr’?’ + 2wr, — w7,), 
(3.12) 
M, = Br’ — or, — w T,) 
and, by substitution in (3.8), the shear forces 
K, = —Bri" + (a + 2B)wr?’ + (2a + B)w’r, — ew, . 
3.13) 
K, = —arl!’ — (2a + Bart’ + (a + 28)wr, + Bw, . 
Introducing them in (3.7), we have 
Bri” — 2a + Burs!’ — 22a + B)wr?’ + Wa + Bw? + Bur, — k, = 0, 
(3.14) 


ar,” + 2a + B)wrl’’ — Wat 28)wr?’ — 2a + Bw, + awr, — k, = 0, 


i.e. a system of 8th order for r, , 7, . 
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Since the rod is loaded only by its inertia forces, 
a = — por. ’ k, = — pi, ° (3.15) 


Substituting (3.15) in (3.14), we obtain the partial differential equations of the flexural 


vibrations. Taking 


r, = U cos xi, r, = Vcosxt (3.16) 
they yield, for the eigenfunctions, the differential equations ‘ 
BU’ — Aa + B)wV’"’ — 2(2a + B)w’U” + 2a + Bw V’ + (Bw* — ux’*)U = 0, 
(3.17) 
aV’’ + 2a + Bw’ — 2a + 28)w°V"" — Wa + B)w*U’ + (aw — ux’)V = 0. 
If the coordinate z is replaced by the angle of rotation of the corresponding section, 
g = w, (3.18) 
we have finally, denoting derivatives with respect to ¢ by primes, 
BU’ — 2(a + 8)V’"" — 2(2a + B)U” + Aa+ pV’ + a( 1 - wu = 0, 


(3.19) 
aV'’ + 2Aa+ 8)U'”’ — 2a + 28)V" — AZa+ BU’ + al 1 -~ uly = 0. 
\ aw 


4. Boundary conditions. The eigenfunctions are determined by the differential 


equations (3.19) in connection with the boundary conditions. In the case of a rod built 


in at the lower end the displacements r, , r, of section g = 0 are zero; thus, by (3.16), 
U= V=0. (oe = 0). (4.1) 
Besides, for ¢ = 0 the tangent of the deflection curve is vertical; hence, by (3.11), (3.16) 
and (4.1), 
U’' = V’ = 0, (¢ = 0) (4.2) 
primes indicating derivatives with respect to ¢. 
If 1 denotes the length of the rod and @ = al its total twist, the bending moments 
(3.12) and the shear forces (3.13) vanish for ¢ = ¢; hence, 


g(U" — 2V’ — U) =a(V" +2U’-—V)=0, &=¢) (4.3) 


and 
B(U'"’ — 2V" — U’) =a(V’”’ + 2U” — V’) = 0. (oe = ¢). (4.4) 
The eigenvalue problem defined by (3.19) and (4.1) to (4.4) is somewhat similar to 
the one of the corresponding buckling problem’; yet, since its order is twice as high, 
it is more complicated. 


5. Solution. The fundamental solutions of (3.19) are of the type 


U = A exp (Ag), V = B exp (Avg). (5.1) 


7H. Ziegler, Schweiz. Bauzeit. 66, 463 (1948). 
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Substituting them in (3.19), we obtain 


ax — 2(2a + B)rX* + a(1 -_= ) | — 2a@+ BAX’ — DB = 0, 


3 a 
(5.2) 
Wy - RY? ae L 92))2 BK _ 
Zia + B)AA" — IDA +] ar Zia + 28)X° + all — 2)1b = Q 
\ aw 
and the characteristic equation 
; 1, 1\ pr’ 
r* + 4° + [ - ( 4 ) a hs 
\Q B ®@ 
| | 2 : "es 2 a 
|. 2| 2 }- 3(2 +. ) a iS re (1 - m1 - us’) = 0). (5.3) 
a B/ @ aw Bw 
By (5.2), any one of its eight roots \, yields a ratio 
B, rs — 21 + 2(a/B))AzZ + 1 — (ux? /Bw') 54) 
= 2S . (9d. 
A, 2|1 + (@/B)JAcQAR — 1) 
and therefore one of eight fundamental solutions. The general solution is 
8 8 
U = > A, exp (A,9), V = DB, exp (rw), (5.5) 
k k=1 


the coefficients A, , B, being determined by (5.4) and by the boundary conditions 


(4.1) to (4.4), i.e. by 


> 4A. =0, B > [Ai — 1A; — 2,B,] exp(d) = 0, 
1 


1 


DA, = 0, BD [Ai — 1) Ax — 2B, JA, exp (Ad) = 0, 
1 1 


(5.6) 
>~B. =0, a@ > [zi — DB, + 2,A,] exp (Aud) = 0, 
l 1 
> rA.B. = 0, a >> [Ak — IB, + 2A.) exp (Ad) = 0. 
1 1 
In order to simplify the problem, we restrict ourselves, by the limiting process 


8 — o, to the most interesting case of a perfectly narrow cross section. If this process is 


carefully applied, (5.3) reduces to 


(A? + 1)* — (A* — 60? + 1) “4 = 0; (5.7) 
aw 
(5.4) yields 
a 3 a 
B, ae A, , (9.8) 





LATERAL VIBRATIONS OF TWISTED RODS 169 


1954] 


and the system (5.6) becomes 


& A, = 0, y 2 Os tM, totes ae) exp (Ax@) A, = 0, 


8 8 2 2 2/ 4 
ps AA, = (), + Os x ©. tae oe ) Ai exp (Ap) A, == 0, 
1 hs 
(5.9) 
> -+— i, ae me = exp (Ao) Ay = 0, 
\ 1 Ai 
> (2 — 1)A, = 0, > (Ai + 1)? exp (A.¢) A, = 0. 


1 


Equating the determinant of (5.9) to zero, we observe that (5.7) is solved by pairs of 
roots, \; = —A , \y = —As. We further make use of (5.7) and of the fact that 
any row may be supplemented by a linear combination of other rows. Thus, we obtain 


the eigenvalue equation 


l 0 
0 Ay 
0 L 
Ay 
IN 0 
1 + A;)* cosh Ao (1 + Xj)’ sinh Ayo = (, (5.10) 
1+)’. ( ji 
sinh \,\¢ a cosh do 
Ay 
3—-XN i--.. 
+2! cosh A\@ ; + 2 Sinh dio 
1—3y_. 1-3; 
Xu. + sinh \,¢ i, + 22) cosh A, 


the three remaining double columns being obtained from the first one by substituting 
Nes Das. te Or Re 
. ~— - ° ° . 
In order,to evaluate (5.7) and (5.10), we introduce, in place of the natural frequency 
x, the quantity 
ux’ l' 


4 
nm = , 
a 


(5.11) 


denoted by k‘/* by some authors. According to 8. Timoshenko’, the eigenvalue equation, 
. ’ . . . 
in the absence of twist, takes the simple form 


1 + cos m* cosh m* = 0, (5.12) 


8S. Timoshenko, Vibration Problems in Engineering, Van Nostrand, New York 1935, p. 234. 
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yielding the eigenvalues 


m* = 1.875, m>$ = 4.694, m= = 7.855, m= = 10.996, --- (5.13) 
Owing to (5.10) and the relation @ = wl, the constant in (5.7) is 
uK m* 
4 ~ 
j=-i=7D; (5.14) 
aw i) 


hence, (5 ean be written 


(1 + A*)* — (1 — 6d’ + d*)p* = O. 5.15) 


In order to obtain the natural frequencies 


m? (a\'” wa 
= i (5.16) 
C \u 


in the presence of twist, it is best to start from a given value of p and to calculate the 
Introducing them in (5.10) and solving for 


corresponding roots d, , «++ , Ay of (5.15). 
, yielding 


¢, we obtain a sequence of values of the total twist, every one of them, by (5.14 
a value 
m = pd. (5.17) 


By repetition of this process (cumbersome since the real or complex character of the 
roots \, depends on the magnitude of p), the functions m,(@), m.(@), --- can be evaluated. 

6. Computer work. The process described above has been carried through by means 
of the sequence controlled computer mentioned in 1. For any of 12 values of p (chosen 


in the interval 0.06 < p* < 26,000) 9 to 18 values of the determinant (5.10) were com- 


puted. The roots \, of (5.15) were introduced right from the start; the angular functions 
of ¢ near the probable zeros were taken from Tables and transferred to punched strips 


for use by the computer. The determinant was evaluated by the method of Gauss and 


Banachiewicz . Controlled by a cyclic program, the computer worked nearly without 
supervision, printing the results from which the zeros were obtained by interpolation. 


Tables 1 to 4 give the values calculated for m,(¢), --- , ms(@), together with their 


TABLE 1 


Fer) | m() m?(¢) 
0 1,875 3,5156 
0, 1477 1,875; 3,51s 
0,6003 1,879 3,53; 
1 ,6209 1,903>5 3,62; 
2, 4328 1,934; 3,743 
2,9273 1,957¢ 3, 83. 
1 0638 2,011. 4,045 





See R. Zurmiihl, Matrizen, Springer, Berlin, Géttingen, Heidelberg 1950, p. 249. 
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TABLE 2 
° m(¢) m3(¢) 
0 4,694 22,03 
0,3599 4,570 20,8» 
0,6775 4,3l1s 18,64 
0, 9566 4,062 16,50 
1 , 2230 3,823 14,6; 
1 ,5869 3,545 12,59 
2,1245 3,223 10,35 
2,5730 3,02, 9,1. 
2,7629 2,950 8,70 
3 ,0656 2,85; 8,14 
3, 4585 2,750 7,56 
3,9630 2,650 7,02 
5.0854 2,517 6,34 
TABLE 3 
| 
¢ ms(¢) m3 () 
0 7,855 | 61,7 
0,6019 | 7,64, | 58,4; 
1,1400 7, 266 52,79 
1 ,6225 6,88, | 47 , 4¢ 
2,0898 6,54, 42,75 
2, 7348 6,115 37 , 40 
3,7083 5,62; 31,66 
4,5604 5,35. 28 , 67 
4.9419 5,277 27 8; 
5,5909 | 5,20; 27 ,07 
6, 5430 5,20, 27,05 
7,9359 5,30, 28,17 
10,857 5,373 28 , 87 
| 
TABLE 4 
od m(¢) m2(¢) 
0 10,996 120,91 
1,598 10,185 103,73 
2,265 9,61; 92,45 
2,890 9 04s 81,8: 
3,714 | 8,305 68,97 
4,873 7,394 54,67 
5,805 6,816 46 , 46 
6, 1975 6,615 43,79 
6, 8314 | 6,357 40 , 42 
7,7011 6,12, 37 , 50 
8,981 6,00. 36 ,0; 
11,888 | 5,88, 34, 62 
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squares which, according to (5.16), are proportional to the 4 lowest natural frequencies. 
Fig. 3 shows the corresponding m(¢)-curves, points taken from Tables 1 to 4 being 


marked by small circles. 
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Fic. 3. Frequency curves. 
Circular Frequencies «; = (m?/l?)(a/p)”? 
1: length of the rod a: flexural rigidity : mass per unit length : total twist 
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For higher modes, the numerical work is complicated by the appearance of differences 
of large numbers. Yet the computer (designed for 6 significant figures) would yield 


more than the 4 curves of Fig. 3. 
7. Expansions. For large values of p, i.e. for small values of the total twist ¢, the 


roots of (5.15) can be developed in the series 


— 5 : - 
A, = (o+2---), dl ee a Ma = V2F1-:-. (7.1) 


The corresponding expansion of (5.10) is 


(1 + cos m cosh m) 


II 


3 0. (7.2) 


p 2 


I , 13 ; ; 
+ -3| 9sin msinh m — = m(cos msinh m + sin m cosh m) | --- 


Developing also 


] ie 
m= m* + > m** «+s , (7.3) 
) 
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we obtain 


* * * 
m m m 
(0 — 13 = cot m=) cot “9? 


~ 


m** = (7.4) 


* * * 
-(9 + 13 “ tan me) tan 
according as m* is an odd or an even root of (5.12). 

Evaluating (7.4) and (7.3) for a number of values of p (chosen sufficiently large) 
and making use of (5.17), we obtain approximations of the frequency curves m,(¢) for 
small values of ¢ (marked, in Fig. 3, by small squares). For ¢ = 0, the functions m, are 
stationary; with increasing ¢, however, m, increases slightly while m, , m; , --- decrease. 

For small values of p, the roots of (5.15) are 





° ? Dp 7 2 
Ae = Ul Fq-:::), As4 = Fit Qe--:. (0 = ba = 725). (7.5 


The corresponding expansion of (5.10) is 


> q ° ° 2 = 
(1 + cos n cosh n) — ; (cos n sinh n — sin n cosh n) = 0, (7.6) 
where 
_— ne oe 
n= wa (4.4 
Developing 
n=n* + n**q---, (7.8) 
we find that n* satisfies (5.12) and that n** = + 1/2. Thus, (7.7) and (7.8) yield 
9 
= 7 m® «os, (7.9) 
6 F 1 


It follows from (7.9) that, for ¢ — ©, two frequencies at a time tend to the same 
limit. The corresponding values of m, according to (7.7), are 


M, = VW2n* (7.10) 


or, evaluated by means of (5.13), 


M, = 2.230, M, = 5.582, M; = 9.341, ---. (7.11) 


According to (5.16), these limits correspond to the natural frequencies of an untwisted 
rod of flexural rigidity 2a. 

In Fig. 3 the expansions for ¢ > @ (like those for ¢ — 0) are marked by small squares. 
For the fundamental mode, the expansions for ¢ — 0 and ¢ — © approach each other 
closely. Evidently, they would suffice for many practical purposes. As to the higher 
modes, however, there remains a gap of increasing width between the two expansions 
(filled by the computing process described in 6). Particularly the expansion for ¢ — 0 


is reliable only for very small values of the total twist. 
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ON THE TIMOSHENKO THEORY OF TRANSVERSE BEAM VIBRATIONS* 


BY 
C. L. DOLPH 
University of Michigan 


1. Introduction. The classical one-dimensional theory of flexural motions of elastic 
bars is based upon the Bernoulli-Euler equation 


7 oO y* a°y* 
EI ax" + p Fy 0, (1) 
where 1*(z, t) is the displacement at time ¢ of a point on the central axis. Equation (1) 


has long been known to be inadequate for the higher modes of vibration. Although 
tayleigh [1] introduced the effect of rotary inertia, the one-dimensional theory did not 
adequately explain the higher modes until Timoshenko [2], [3], extended it to include 
the effects of transverse shear deformation. The Timoshenko theory is based upon the 


following equations: 


ie,  _ 2 
or + ox 0, - 
” * 22% 
om = + Io £z.. + sia - 0, (3) 
Ox ot 
op* M* | 
ox ~—«C«&EW’ ™ 
V* = —KAGao*, (5) 
90)* 
= = at + B*, (6) 


In these equations y*(x, t) represents the transverse deflection, 1/*(x, t) the bending 
moment, V*(z, /) the transverse shear force, 6*(x, t) the neutral axis slope due to bending, 
and a*(x, t) the neutral axis slope due to shear. Also, EJ represents the flexural rigidity, 
KAG the shear rigidity, and J, the so-called rotary moment of inertia. 7) is equal to 
ol /A where p is the mass per unit length, J the area moment of inertia, and A the cross- 
section area. The first two of the above equations are stress equations of motion, the 
second two are stress-strain relations, and the last is a geometric strain-displacement 
relation. 

While these equations may be derived in a plausible fashion from the laws of dynamics 
and of elementary strength of materials, their ultimate justification, as well as the range 
of their validity, must be deduced from the general equations of elasticity. This has 
recently been accomplished by Mindlin and Hermann not only for beams but for many 
other approximate equations of elasticity. While these results have not yet been pub- 
lished, they were presented at the Maryland Symposium in Elasticity in the spring of 
1952. Some hints as to their method may be gleaned, however, by reference to Mindlin 


(4]. 


*Received March 21, 1953. 
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The present paper presents a derivation of some of the consequences of the 
Timoshenko theory. These same results could be deduced by use of the basic approach 
of Mindlin and Hermann. However, since the Timoshenko theory has formed the starting 
point of many investigations (e.g., [5], [6], [7], [8], [9]) an elementary approach based 
directly and solely on the equations of that theory will be employed here. 

The simplicity of the treatment presented here results from the fact that, apart 
from the last section, it is based directly upon the equations which result from the above 
when a* alone is eliminated. This is in contrast to the more usual treatments which are 
based upon the involved displacement equation obtained by eliminating all variables 
rom the above set. Some idea of the complexity introduced by this latter 


# 


except y* (x, 
procedure can 
beams which is known to be the following: 


EIp ) a°y* I,p d'y* o°y* 
es ow oe oe oO eee igo me’), ( 
(7 , KAG/ dx dt KAG at ? 


: 
i 
be gathered from an inspection of the displacement equation for uniform 


_ 0 y* 
EI = 


ax" 


~I 
— 


ot 


The main consequence of the approach in this paper is the fact that the discussion of 
orthogonality, initial value problems, problems involving external forces and moments, 
and problems involving time dependent boundary conditions for the Timoshenko theory 
can be mad 
ry. In particular, for the uniform hinged-hinged beam, a complete analytical 
n be given (much of it will be presented here) to all of the above problems. 


le to parallel step by step the familiar treatment of the classical Bernoulli- 


1 
} 
il 


Euler the 
solution ca 
Moreover, for non-uniform beams, the computation of characteristic values by analogue 
simulation or by other automatic means is much more direct if the system approach is 
used and the results of this paper can be applied immediately (cf. [9]). 

There are of course important points of difference between the Timoshenko theory 
sed on equation (1), or upon theories which neglect either the correction 


and that b 
The most important of these is the existence in 


due to shear or rotary inertia. 
the Timoshenko theory of two sinusoidal modes of different frequencies corresponding 
to the same spatial factor. These two sets of modes are, moreover, orthogonal to each 
other in terms of the two term orthogonality relation which forms the kernel of our 
method. The existence of the second set of modes also makes it possible to satisfy the 
four initial conditions required by the Timoshenko theory. 

Finally, the existence of the second set of modes makes it possible to demonstrate the 
completeness of the set of particular solutions, sinusoidal in the time, for the above 
named problems of the Timoshenko theory. 

2. Separation of variables and the orthogonality relation. The separation of variables 
process can be applied directly to (1) but not to either (2-6) or (7). Thus, for example, (7) 
will separate under the assumption that y*(x, 1) = y(x) T(t) if and only if one of the 
quantities 7, , 1/KAG vanishes. All will, however, admit solutions sinusoidal in time. 
Thus if we set 3/*(z, 2) (x)e'**, V*¥(x, t) = V(x)e'**, M*(2, t) = M(x)e'*’, B*(z, t) = 
B(r)e’*' and eliminate a* from the system (2-6), it reduces to the following: 


dV 2 (OQ 
- = poy, (3) 
da 
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dB a 

ai” a” (10) 
dy _ - ae 
dx °° + kag = ° (11) 


In order to discuss the orthogonality of solutions sinusoidal in time for a beam of length 
L we note that integration by parts furnishes the following Lagrange identity appropriate 
for the system (8-11) in the usual way: 


oat dV, (dM, a (# M,) 
[iL SB - ol — v,) + aol - 3 


(2 — 9. Mee] — fin M2 — (la — v,) + an (8 — 3) 
Abe Bi 4 as | E dx B dx Vs) + M, dx EI 


_fdyp vs : . - : 
F : — | 1. : ‘= x — fo ms —_ 9 . 2 
J (2 2 + aa | dx ve 1 — BM, + MB, J a (12) 
Thus if y, , V; , 8; and w, , (¢ = 1, 2), represent two solutions of (8-11) with boundary 


conditions that make the right-hand side of (12) vanish, the following orthogonality 
relation is obtained, independent of whether the beam is uniform or not: 


L¢ 


~L 

(oi — 2) | Loyrys + ToBsB2] dx = 0. (13) 
Although many boundary conditions are possible, this discussion will be limited to the 
consideration of free ends where M* = V* = 0, hinged or simply supported ends where 
y* M* = 0, and built-in ends where y* = 8* = 0 (although some authors mistakenly 
use y* = dy*/dx = 0 in the Timoshenko theory)’. Thus, from (12) it is apparent that 
(13) will hold for: (a) the free-free beam; (b) the cantilever beam (i.e., one end built-in, 
the other free); (c) the hinged-hinged beam; (d) the built-in, built-in beam; and (e) the 
built-in, hinged beam. 

If J, is zero, and w; ¥ w , the relation (13) reduces to 


~L 

| pyiy. dx = 0 (14) 

Jo 
which is identical to that of the usual textbook theory. This shows that the displacements 
associated with different. modes are orthogonal with respect to the weighting function 
p\ DP 
Even if J, is not zero, the relation (13) is readily interpreted. As we shall see in §6 
the solution of the initial and boundary value problem will, by our method, reduce to 
the determination of a set of constants C, (n = 1, 2, ---) such that two given functions 
f(x) and g(x) admit the simultaneous expansion 


f(z) = > C,Yn(X), g(x) = > C,B,(x). 


1For a discussion of this point as well as of the orthogonality relation for the cases where either the 


shear or rotary inertia term is neglected, see [10]. 
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In order to show that this is possible in virtue of the relation (13), we introduce two 
mutually perpendicular unit vectors e, and e, and consider formal vectors of the form 


A(x) = a,(x)e, + a.(x)ez , 


with a scalar product defined by 
aL . 
| A(x) - B(x) dx = | [pa,(x)b,(x) + Ioay(x)b.(x)] dx. (15) 


If in particular we let 
K(x) = f(x)e, + g(x)e. , Z(x) = y,(x)e, + B,(x)e2 , 


then we must determine the constants {C,} in the vector expansion 
K(z) = > C,Z,(2). 
p=1 


Taking the scalar product in the sense defined by (15) we find that 
i J | pf 2) us) + Tog(x)B,(2)] dx Bi ag i Se ich, (16s 
fo [ofy,(x)}° + Io{8,(a)}"] dx 


The situation here is somewhat special in that for a general fourth order system one 
would expect that it would be necessary to introduce four mutually perpendicular unit 
vectors and a scalar product of the form 


aL 


| Z,(x) - Z.(x) dx = | [oyiys + ToBiB2 + AMM, + BV,V,] dz. 


Because the last two terms are absent from (13) one cannot in any real sense say that the 
vector solutions y, , V, , M, and 8, , corresponding to w, , are orthogonal to each other 
for different k’s but only that the first two components of different solutions have this 
property. From this point of view the relation which obtains if J, vanishes is even more 
special in that then only the first components are orthogonal.” 

3. Normal modes of the uniform beam. If the beam is uniform the system (8-11) 


can be treated by standard methods. Thus if y(x) = ye*, M(x) = Moe’, V(x) = 
V,e7, B(x) = Boe” where yo , Vo , My and ) are constants, and if the notation a 1/EI 


and b = 1/KAG is used the system (8-11) reduces, after division by e”*, to 
AVo — pw Yo = O, AM, — Vio + Iw Bo = 0, 
AB, — aM, = 0,7 AYo — Bo + LDV, = O. 
This system can have a non-trivial solution if and only if the determinant of the co- 


efficients vanishes; explicitly, if 


h* + (alo + bp)dA*w? — apw’ + abplw* = 0, (17) 
from which it follows that the roots \ satisfy the equation 
Dn? = —(aly + bp) + {(aly + bp)*w* + 4apw’ — 4abplw'}'” 


2M, V may however be eliminated as in §8. 
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There are several possible cases associated with this characteristic equation. If w = 0 
then all of the roots are zero and the system (8-11) simplifies in that its right-hand 
members vanish. It is readily verified that the solution in this special case is the following: 


( yo 0. 2 03 a 6CY 
7! G + Cox + Cyr" + Ce’; ili — . 
— 4 6C or 0 6 rn 0 0-5 
M = ls (c: — b cr) + 2C3x + 3C%x’. 
a a 
If w = 0 but w» = 1/b/, then \ = 0 is a double root and the other two roots are 


purely imaginary and conjugates of each other. This case cannot be admitted into the 
theory of the system (8-11) since it implies that the relations for V and @ in terms of y 
that can be deduced from (8-11) are meaningless (cf. [9]). 

For reasonable physical beams, the condition w” < 1/b/, is valid for the frequencies 
of interest. Under this assumption we obtain the four roots: 


Ai. Ae = 4; As, Ae = te, 
where 


2 J = {—(al, + bp)w* + [(aIy — bp)’w* + 4apw*)'’*}'”. (18) 
€ 


From this it is readily deduced that the most general solution of the system (8-11) must 
be, for fixed w, of the form 


y C, cosh yx + C, sinh ya + C3 cos ex + Cy sin ez, 


, pPwirn -: ’ ee ’ 
V = _ {C, sinh yx + C, cosh yx} + . {C, sin ex — C, cos ex}, 
i 
(19) 
bow +7 ., ’ bw —eé : 
M = = — {(, cosh yx + C, sinh yx} + == {C, cos ex + C,sin ex}, 
a a 
bow Tr o F (sy ° y ? b pw” — € if’ . y t 
B C, sinh yx + C, cosh yx} + ———— {C; sin ex — C, cos er}. 
€ 


I 


t. Normal modes of the uniform hinged-hinged beam. In order to find the modes of 
oscillation in the cases of physical interest, boundary conditions need to be imposed on 
the equations (19). By far the simplest case and the only one which seems capable of 
complete analytical treatment without approximations is that of the hinged-hinged beam 


where y* = M* = 0 holds at x = 0 and x = L. Imposing this boundary condition at 
x = 0 we obtain the equations 
b ? - b RS _ Pa " 
C.4+C, =&, oe C:.% tin —C, =0. 


Since the determinant of this system has the non-zero value —(y° + é’)/a, it follows 
that C, = C, = 0. It is convenient at this point to simplify the formulas by selecting 
the unit of length so that the beam is of length 7. The formulas for a beam of arbitrary 
length can be obtained by replacing yx and ex by yL and eL respectively in the relations 
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(20, 21) and n by nz/L in all succeeding formulas. The boundary conditions at x = 7 


yield the equations 

C, sinh yr + C, sin er = 0, 

(20) 
(bow + y7°)C. sinh yr + (bow — €)C,sin er = 0. 

In order to have a non-trivial solution the determinant of this system must vanish. Thus 
one obtains the equation 

(y*> + €) sinh yx sin er = 0. (21) 
Since y and «¢ as defined by (18) are real, the only roots of interest are simply those for 
which 
y ~ 0, e=N, ieee OD, 2 Ze =), (22) 


i 
for if y and e were both zero, then w must be zero and the special solution of the system 
(8-11) must be employed. It is readily verified that there exist no non-zero constants 
Cy (k = 1, 2, 3, 4), which permit the boundary conditions to be satisfied for w = 0. On 


k = 9 ay 
the other hand, the combination y = 0 and e = n cannot occur, since this would imply 
that b/w = 1. 

The corresponding spatial parts are therefore of the form 
te [bpw —n\. 
y = C,, sin nz, M= oe , sm nz, 
a 
e ; pw _(n? — bpw 
y =o { — COS NX, B= c, COS NX, 
\ n | n 
forn = 1,2 


It should be noted that for the uniform hinged-hinged beam the assumption that 
bIyw < 1is really unnecessary, for if bJ,w > 1 it follows that all four roots of (17) would 


be imaginary so that the secular determinant would involve the product sin yz sin em. 


This product does not introduce any new spatial eigenfunctions. 


The values of w which correspond to « = n are obtained by setting \ = ze = in in 
equation (17) and solving the resulting equation for w. Thus, from 
‘ | 22 + (al, + ben 2 n' o 
o- a w + -=0 (23) 
abpl, abpI, 


we obtain 
lap + (al, + bp)n*] + {lap + (al, + bp)n*} — 4abpIyn*}' Z 
x 2abpl, 
so that corresponding to the spatial functions (22) there are two distinct values of w’. 
Let the roots be denoted by + w,, and + w., . We shall assume | w,, | < | wy», |. For 


fixed n, we therefore obtain particular solutions to the system (8-11) of the form 
y = sin nx{E;,, cos w,,t + F;, sin w;,t}, 
V = —cos na[—pw;,/nj|[E;, cos w;,t + Fi, sin wnt], 


(24) 
M = sin na[(bpw:, — n’)/al[E jn COS Wint + Fin Sin win], 
8B = —cos na[(bpw;, — n°)/n\[Ein COS wint + Fj, Sin wint]; 


(¢ = 1, 2:n = 1, 2, ---). 
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Consideration of the second mode which occurs for n fixed seems to have been entirely 
overlooked by Timoshenko [2] in his original paper, although he treated only the uniform 
hinged-hinged beam. 

For different values of n it is clear that orthogonality results, since the general relation 
(13) is satisfied trivially inasmuch as the integrals 


aT ar 
p Yio dx, I, | 8,8. dx 
70 “0 
vanish separately. We will now investigate how the relation (13) is satisfied for the same 
n but for the two modes which correspond to w,, and w., . In this case, Equation (13) is 


explicitly 


“0 


| (py Ye + IB,B2) de = K | {o sin’ NX 


 g 992 2, 2 2 
+ =$ [b?p*w?,wen — bon*(win + won) + n*] cos” nz} dx, 


n 
where K is some arbitrary constant. 
The coefficient of cos’nx in this expression can be simplified by noting that it follows 


> 


from (23) that 


ap + (al, + | bp)n? 


Win + Won = 
abplI, , 


and that 


Thus 
To 4,222 2 2 2 2 4 
n° [b P WinWen — bpn (win + Won) +n ] _ me 
so that 
| (pyiYo + Io8,82) dx = Kp | (sin? nx — cos’ nz) dx = —Kp | cos 2nx dx = 0. 
“0 “0 


Thus the modes are orthogonal in accordance with (13) if they correspond either to 
different n’s and the same w or to the same n and to different w’s. It is necessary to use 
both sets of modes for a fixed n in order to obtain a solution to the general initial and 
boundary value problem. 

5. Comparison of universal curves for hinged-hinged and free-free beams. In order 
that the shear and rotary inertia frequency effects may be estimated for any uniform 
hinged-hinged beam, it is convenient to introduce dimensionless parameters. Such a 
set of parameters has already been introduced by Howe et al. [9] and an electronic 
differential analyzer has been used to compute a set of curves showing the effects of 
rotary inertia and shear on a uniform free-free beam. Since these authors have kindly 
granted permission to reproduce their curves, it seemed expedient to use the same param- 
eters and to compute the corresponding values for the uniform hinged-hinged beam 
for which the exact theory presented here is available. 
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The dimensionless parameters introduced by Howe et al. were defined as follows: 


S = EI/KAGL’?, N=KG/E, d= olL"(p/EI)"””’. 


If the effects of shear forces and rotary inertia are both negligible, then S = N = 0. 
This will be the case if the length of the beam is very long compared with its thickness; 
such an idealized beam is often termed “‘infinitely long’’. If A, is used to denote the value 
of the appropriate mode of the “‘infinitely long’ beam, the ratio \/A, is a convenient 
measure of the change in frequency caused by values of S and N different from zero. 
For the uniform hinged-hinged beam, \, may be computed directly from equation 


(23) which may be written in the form 


NSN — [1 + S(1 + N)m*r' Jr, + m'x* = 0, (m = 1, 2, -->), 25) 
so that 
Aue = (2N)'7[S-? + S11 + N)m*r” 

+ {{S-?+ S11 + N)m’n’) — 4S°-?m*a'N}'?]'”. 26) 
The corresponding values of \, are computed from (25) by setting S = N = 0 so that 

Ny = mr’. 
The smallest values of 100 \,/A, for a particular mode are plotted versus | S in 
Figure t for values of NV 0, 0.1, 0.2, and 0.3. The corresponding values for a uniform 
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free-free beam as given in Howe et al. [9] appear in Figure 2. Figure 3 is a similar graph 
for the values of \,/A, as determined by the plus sign in equation (26) for the uniform 
hinged-hinged beam. 

These curves may be used to estimate either separately or simultaneously the effects 
of shear and rotary inertia upon the frequency of vibration. It is apparent from its 
definition that S increases with an increase in flexural rigidity, and decreases with an 

















1954 TIMOSHENKO THEORY OF TRANSVERSE BEAM VIBRATIONS 183 

























































































100 as/aq 
100 = 
ee = Be —> _—<—— 
LT] —— wa“ 
Ces 
90 Cs LA] 
Vint 
| 
80 3 A 7 
— | 
” 
(oe | 
| | | 
70 |}— ! 
| | | | 
al 
4 | LL vs 
10 1000 5000 
Fic. 2 





10 100 1000 5000 
Fie. 3 


increase in shear rigidity or an increase in length. It is therefore a parameter capable of 
describing the relative importance of shear rigidity for beams of given flexural rigidity 
and of given length. On the other hand, the expression for N may be written as 


N = KAGI,/EIp (27) 


and is therefore seen to increase with an increase in shear rigidity or rotary inertia and 
to decrease with an increase in flexural rigidity. Thus, N is a suitable parameter to measure 
the relative effect of rotary inertia for beams of given shear and flexural rigidity. Alter- 
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nately, since J, is directly proportional to the density, it may be used to measure the 
relative importance of the size of the cross-section. 

6. The initial value problem for the hinged-hinged beam. 


An examination of the 
system (2-6) shows that the general initial value problem requires the specification of 
y*(x, 0), y*(a, 0), B*(x, 0), and B%(zx, 0) in addition to the boundary conditions appropriate 
to the type of beam under consideration. Let us write 


y*(z, 0) = f(a), B7( 2. 0) = gi(x), 


(28) 
y*(x, 0) = f(z), B*(x, 0) g2(x). 


9) 


Then, because of the linearity of the system ( 


-6), a superposition of the solution (24) 
will satisfy this sytem and the boundary conditions for the hinged-hinged beam. Such 
a sum is of the following form: 


—— » a sin n2z[E in COS wnt + F;, Sin w; t] 


B*(x, 1) = >> > cos nz(bow:, — n’)n"' 


where M*(z, t), V*(; 





E,, cos w;,t + F;, sin w;,¢] 


|, 
(x, ¢) can be expressed in a similar fashion. In order to satisfy the 
initial conditions (28) 


, the coefficients Ein ‘ F, ~4@ = 1 2n Sy. 


. ,2;n = 1, 2, ---), must be de- 
termined from the following series representation: 


fiz) = DOE, sinne + > E,, sin nz, 


n=1 


g(x) — YE, (bpwin — n’)n™ cosnz — >> E>,(bpw2, — n’)n™' cos nz, 
f(z) = Y nF.) sinnz + D> (2 2,) sin nz, 

n=1 n=1 
J2(2) = pe (winks, )\ bpwin ee 


2) -1 ’ 2 2) - 
n\n cosnx — > (Wonk on)(Dpwon — n')n' cos nz. 


n=1 

Although the coefficients E,,, , E.,, F:, , and F;, can be evaluated in this case directly 
by the usual Fourier methods, the orthogonality relation (13) and its interpretation 
offers the advantage of furnishing these values directly without the necessity of solving 

two sets of algebraic equations. Thus, for example, 
¢ Pe / 4 -1 2 2 
a 2 fo [f:(az) sin nz — Ion 'p (bpwin — n°) gi(x) cos na] dx 
417 oo ; 1 2 q 2 2\2 oe 
T 1+ p n “I,(bpw;, — n’) 
If the usual Fourier coefficients are employed one finds 
9 


7 “a 
Ein = 


fc (f(a) sin nz + n(bpew2, — n°) 'g,(x) cos na] dx 
T 


bp(win — Win) [bp(w2, — n°)’ 
The equality of this with the first expression for E,,, follow directly since 


n/(bpwen — n°) = —I,bpwin — n’)/np, 
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and 


> Dy 


(Down — Win) /(bpwen — n°) = 1 + (n? — dpwin)/(Dpwon — n’) 
= 1+ I,(bpw;, — n’)?/pn’, 


in virtue of the sum and product of the roots w,3 and w3, of (23). 

The rigorous justification of the above expansions as solutions of the initial value 
problem would, in general, involve a proof of the completeness of the set of eigenfunctions 
and a study of the convergence properties of these expansions. Here, because of the 
identity of these expansions with those of usual Fourier sine and cosine series, it is well 
known that the set of functions is complete and many conditions sufficient to guarantee 
the convergence of such series are, of course, known. Moreover, since the general solution 
of the system (8-11) for # = 0 does not admit a special solution satisfying the boundary 
conditions and allowing V* and 8* to be constant, the identity of these expansions with 
those of Fourier series can be used to deduce the fact that the functions g,(x) and g,(x) 
cannot be specified arbitrarily but must satisfy the additional conditions that 


| g(x) dx = 0 for7 = 1, 2. 


“0 


This conclusion is however a function of the boundary conditions. For example, for a 
free-free beam, the solution 8 = constant is permitted by the system (8-11) if w = 0. 
For the uniform hinged-hinged beam it is clear that these constraints on g;(x) are merely 


symmetry requirements since 8 is the slope due to bending. 


7. Beam problems involving external forces and moments. If an external force 
F*(x, 1) per unit length and an external moment G*(z, t) per unit length is allowed to 
act on a non-uniform beam, only the first two differential equations of the system (2-6) 
are modified and they become 

oM* 


i ae" - a°B* 
y T 4 = F*(z, t), —_ a | of 


— '? = G*(z, i). 
ot Ox Ox or” we Gz, 0 


p 

The solution of any problem of this type may be obtained by the superposition of 

solutions of the corresponding problem in which the external forces and moments are 
zero. It is merely necessary to seek solutions of the form 


at at 


of? =x | y'(x, t — a) do, V* = | V'(a, t — a) do, 
(29) 
Mt= | M(z,t—o)dc, pt=| Be, t— 0) do, 
“0 “0 


and to impose the following initial conditions on the functions y’(x, r) and 8’(x, r) where 


y'(z, 0) = 0, 
y’(x, 0) = yi(x, t — o) |--, = F*(z, 0), (30) 


B’(x, 0) = 0, Bi(a, 0) = B(x, t — o) lene = G*(z, 0). 
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It now follows that the functions y’, M’, V’ and 8’ satisfy the free system (2-6) and that 
if they satisfy a given set of boundary conditions then the function y*, M*, V* and * 
defined by (29) will satisfy the same boundary conditions. In order to determine the 
functions 7’, M’, V’ and £’ it is merely necessary to solve a special case of the initial 
value problem corresponding to the given set of boundary conditions and the initial 
conditions (30) as a function of the parameter. 

8. Beam vibrations with time dependent boundary conditions. In order to reduce 
a beam problem for the system (2-6) which involves time dependent boundary conditions 
to manageable form, it seems simplest to eliminate M* and V* from the system (2-6) 


to obtain the equations 


9° y* 4) ay* 
L(y*, 6%) = p 4 + — E AG(s* = au") = @ 
ot Ox Ox 


L,(y*, B*) = — = (zr 26°) + I, oe + K. 1G(a° - ay") = 0 
and the following expressions for M* and V*: 
ut = EI, y* = KAG(s* - oY) 
OX Ox 
Now if y* and 8* are decomposed to the relations 
y*=yst+yr,  B* = Bt + BF 

these become 

L,(y* , 8%) = —Li(y% , 8%), LAy% , 6%) = —L.(yt , 84), (31) 

M* = M* + Mt, V* = V3+ Vt. (32) 


If, following Mindlin and Goodman [11], one sets 
y* = f,(dhi(x~) + f.(ih.(z), B¥ = gi(Ok,(x) + g.(bk.(2), 


where the functions f;(x), g;(x), h;(t) and k;(t) are arbitrary for 7 = 1, 2, it is always 
possible to choose these cight | functions in such a way as to iui the boundary value 
* , 8% , V4, and M¥% to one involving time independent boundary con- 


problem for z 
ditions, known initial conditions, and known external forces and moments as given by 


the right-hand sides of (31, 32). This problem can in turn be decomposed into two 
problems of the types created in §6 and §7 respectively and thus a general solution 
can be obtained by superposition. 
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—NOTES— 


NUMERICAL QUADRATURE OF SOME IMPROPER INTEGRALS* 
By H. SERBIN (Purdue University) 


1. Introduction. Problems in potential theory, such as occur in aerodynamics, may 
be formulated in terms of integrals which are usually singular. An explicit solution of 
such equations is rarely feasible and one is compelled to turn to numerical (approximate) 
methods. One such method, exemplified by Multhopp’s solution [1] of the equation of 
the lifting line, depends on a quadrature formula of the singular integral. In the following 
analysis, we derive several formulae related to Multhopp’s. 

2. Parseval’s identity. Let @ be a variable ranging on the interval (—7, 7), let 
6, ,¢%, = rx/N where r = 0, +1, --- , +(N — 1), N and N isa positive integer, and 


consider the formal Fourier series on the interval (—7, 7) 
f(@) = >. a, cosn@ + Zz b,, sin n@, g(@) = 2. ay cosné + 2. b sin né. 


The N-th section of a Fourier series will be defined by writing 


Vv 


V 1 
fr(@) = > a, cos né + % b, sin né. 
0 1 


Introduce the notation 


[f(0), g(0)] = (2r)~’ | f(Ag(@ dé, [f(0), g(@)]w = (2N)* > f(8,)9(0 


where the summation is over all r. Parseval’s relation is 
Bs ,,le er ’ 
[f(0@), g(@)] = aaj + = z. (a,a, + b,b%) (1) 
-~ n=l 


and the corresponding relation for finite harmonic series is 
N-1 


; ] 
(fv , 9n)w = Qoag + a z (a,a, + b,b/) + ayay . (2) 


aad l 
Comparing Eqs. (1) and (2), one obtains the 


THEOREM 
(fv, gv) = (fv, 9n)w — 4yax . (3) 
The theorem establishes a means of replacing an integral by a finite sum. When 
fw and gy are replaced by integrable functions, then the right side of (3) is an approxi- 
mation to (f, g). 
For a fixed N, the functions fy(@) form a linear manifold in which two functions 
fw(@) and gy(@) are equal if and only if for every r 


fv(9,) = gy(6,). (4) 


*Received July 15, 1953. 








1954] H. SERBIN 189 


Equality will be used only in this sense below. We shall be concerned with linear trans- 
formations 7’ which transform this manifold into itself. 


As an example, consider the identity transformation 7 defined by 
T[fw(0)] = fr(O). 
Since T(1) = 1, T (cos n@) = cos n6, T (sin n@) = sin né@, one can verify, using the 


Fourier orthogonality relations, that 7’ may be represented by writing 


Tfx(0) = [ +2 > (cos n@ cos np + sin né sin nd), sald) 


N N-1 
[Thv(0) lone, = E +2 2 cos n6, cos nd + 2 p sin n6, sin ng, ful@) 
1 1 
Applying the theorem above, 


V N-1 
[Th v(@)]o-0, = [ +2 > cosné, cosnd + 2 > sin n8, sin ng, ful) | — (-l)’ay. 
1 1 


N 


3. An integral in the theory of airfoil sections. An integral which occurs repeatedly 
in the theory of airfoil sections defines a transformation 


° 


TIs(@)] = 2m)" [| cot 5 (6 — g)F(@) ae (5) 


where the integral is to be understood in the sense of “principal value’’: 


a® nb—e ® 
iN = Lim i + rt 


e—0 


If f(@) is regarded as being defined on the unit circle then T[f(@)] is that harmonic 
conjugate of f(@) which has average value zero. One has 


T(cos n@) = sin n6, T(sin né) = —cosné. 


Hence one has the following representation of 7 
N 
T(fv(0)] = | 2 > (sin n@ cos n@ — cos né sin ng), su(@) |. 
1 


But [sin N¢, fy(¢)] = 0 and the summation of the second term under the bracket need 
extend only to N — 1. Applying the theorem and noting that sin N¢@, = 0, one obtains 


N 
Tl fx(@)] = 2 a, (sin n@ cos n@ — cos né sin nd), Jul) | 
N 


’ 
= 2 >» sin n(@ — @), fu(@) | 
1 N 
Use the identity 


N 


> sin n(@ — ¢) = cot s (6 — ¢), 0 


according as 6 — ¢ is an odd or even multiple of r/N. Then 
Tfx(0)] = [2 cot 3(6 — 4), fr@)]x , 
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where the prime indicates that the summation on ¢ is taken over those values of ¢ which 
differ from @ by an odd multiple of r/N. Explicitly, 
T(fxv(0)Jo-0, = N~ >>’ cot 3(0, — ¢,)fx@,). (6) 
r odd 
4. Example. To illustrate the preceding result, consider an infinitely thin stationary 
airfoil y = y(x), —1 S x S 1, immersed in a perfect fluid of uniform remote velocity 
(1, 0). Let a(x) represent the function dy/dzx and let y(x) be the vorticity per unit chord. 
Then one has the equation [2, p. 88] 
l 
a(x) = (2r)7’ | (x — 2’)""y(x’) dx’. 
J—] 
Make the transformations x = cos 6 , 2’= cos @ and multiply through by sin @. 
Using the same functional symbols a and y, one has 


a(@) sin 6 = —(2mr)7’ | sin 6(cos 6 — cos 6’)"'[y(6’) sin 6’] dé’ 


“0 


; i ° l : en ‘ He 
= (49r) ~ | | cot == 6") + cot 5 (6 + v) [Ino sin 0”) a0" (7) 


Extend the range of definition of a(@) to the interval (—7, 0) by defining a(@) as an even 
function and extend y(@) as an odd function. Then (7) becomes 


a7 


2a(6) sin 6 = (2x) | | cot : (9 — 0) |ivo” sin 6’] do’ = T(f(@)] (8) 
where f(@) = y(@) sin @ and 7 has the meaning of Eq. (5). 
In airfoil theory, one invokes the Kutta condition (y is finite at the trailing edge 
x = 1,i.e., @ = 0) so that f(0) = 0. Applying Eq. (6) and reversing the argument leading 
to (7), 
TIf(0) lene, = N™ py cot (6, — ¢.)f(@.) 
{ 


(N-1 ° 
Bae at : _— sin @ ; ] ] - 
—=T(f(®]e-0, = N" mao i = ev tan = 6,f(w) ? (9) 
."* aan 6, — cos ¢,° o) +5 ® lat Pn) | 
where ¢., = 0 or 1 according as N — ris even or odd. Comparing (8) and (9) and simpli- 
fying, one obtains an approximate equation 
k= ! 
—a(@,) = N™* py (cos 86, — cos ¢.) f(¢,) = = e,v(1 + cos 6,) Sox), 6. + 0, =. 
s=1 “a 
Returning to the original z-coordinate system 
(N-1 1 ) 
a(x.) & —N~") D0’ (a, — 2.) f(z.) — 5 ev(1 + 2,) ‘I De, (10) 


1 ~ 


where | cos 6, | ¥ 1. 

Equation (10) furnishes N — 1 equations in the N unknowns f(z,), --- , f(rw). 
One additional equation follows upon the observation that the left side of Eq. (8) con- 
tains sine terms only up to order N — 1. Hence f(@) must contain cosine terms only 
to order N — 1. Therefore 


ay = (2N)"' >> cos N@,f(6,) = 0, 
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which leads to the equivalent equation 
— f(a) + f (x2) = f(x) teert (3)f(-—1) = 0. (11) 


Equations (10 and (11) furnish a system of N equations in N unknowns. 
As a numerical example, consider the straight line airfoil inclined at a small incidence 


a to the airstream. Then a(x) = —a and Eqs. (10) and (11), for N = 3, give the system 
fe = 3a, —f; + fs = 3a, —fi + fe = fs = 0, 
where f, = f(z,). Solution gives f, = a, fz = 3a, f; = 4a. Since fy = 0, one may inter- 


polate harmonically with respect to @ using N = 3, obtaining 

f(0@) = 2a(1 — cos 8), 

7(0@) = 2a(1 — cos 6)/sin 6 = 2a(1 — 2)'7(1 + 2)”, 
which is the exact result. 

The preceding example is not intended as a substitute for the more direct method 
usually used in solving the problem of thin airfoil sections. However, it does illustrate 
the utility of the quadrature formula (6) in the evaluation of integrals (5) which appear 
in other connections (thick airfoil theory). 

5. An integral in the theory of the lifting line. As the next case, consider the trans- 
formation 


aT 


a6 : ms 
T(f(@)] = (2m) | cot 5 (8 — o)f'@) dd, (12) 
where the prime denotes differentiation. One has, as before, 


™ 
Tf v(0)] = 2 > sin n(6 — ¢), sie) |. 
[ntegration by parts gives 

Tl fxv(0)j = |2 » Bs n cos n(@ — ¢), 56) | 


and application of the theorem 


Tf x(O)Jo-0, = |2 yn cos n(6, — ¢), fu) | — N(cos N@,)ay . (13) 
Use the relation 
ay = (2N)" a (—1)"fx(6.) 
and the identity 
2 > ncosnz = N, —N — ese’ 2/2, N+N’, (x = kr/N), 


according as k is even and different from zero, k is odd, or k = 0. Then (13) takes the 


form, after simplification, 


Tf x(8)Jon0, = 4Nfx(0,) — (2N)~* QU’ esc” $(8, — .)fx(@.). (14) 
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It is interesting to compare (14) with the result of integrating (12) formally by 
parts, disregarding the singularity. The integration leads to a divergent integral 

mir \-1 2 1 ‘ . ~ 

T{f(6)] = —(4m)"' | ese? 5 (0 — o)f@) do. (15) 


v-—f 


The sum in Eq. (14) is a special kind of ‘Riemann sum” corresponding to the integral (15). 
Equation (14) may be applied to derive Multhopp’s quadrature formula. Consider 
a finite wing of unit semi-span immersed in a uniform (perfect) stream of unit velocity. 
Let y and y’ represent the spanwise coordinate; let Aa(y) and ['(y) represent the induced 
angle of attack and circulation, respectively, at station y. 
The induced angle is given by [2, p. 135] 


“ 1 
:' IT /dy’ : 
Aa(y) = (4m) | at A if dy’, (16) 
J-1 Y—- Y 
where 
T(—1) = Pl) = 0. (17) 
Put y = —cos ¢, 7’ = —cos ¢’. Use the same functional notation and extend Aa(¢) 


and I'(¢) into the range —z < ¢, ¢’ < 0 as even and odd functions respectively. In a 
way similar to that used in deriving Eq. (8), one obtains 


4 Aa(¢) sin ¢@ = T|T(@)], (18) 


where T is used in the sense of Eq. (12). Using (17), one can write the sum in (14). 


> By ese +(0, — ¢,)I'(@¢,) = by { ese? $(0, — o,)T@,) + ese” +(0, + ¢,)I'(—¢,)}. 
oft s 
Replace T(—¢,) by —I(@,). Passing to trigonometric functions of double argument, 


one obtains 


py ese” (0, ~~ ,) T'(@,) = 4 > Dy sin 6, sin @,( COs #6. — cos ,) T(¢,). 


— 
Substituting this result in (14) and the result so obtained in (18), one obtains, after 
simplification, 


Aa(@,) = NI(6,)/(8 sin 6.) — (2N)"’ >-’ sin ¢,(cos 6, — cos ¢,) °T'(¢.). (19) 


0«< <N 


Equation (19) is equivalent to Multhopp’s quadrature formula. For numerical applica- 
tions, the reader is referred to [1]. 
6. An integral over an infinite range. Define a-transformation T by 


T(f(x)] = 2" | (x — y)'f(y) dy. (20) 
Consider the related transformation 
T [f(x)] = 3 | (x — y)'f(y) dy, (21) 


where L is a positive quantity. 
Make the transformation 
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and write 

T (f(z) = F(), = fly) = fi). 
Then 

F() =x" | (6 — 8) "f.d) ae. (23) 

Write 
6 — ¢ = 3 cot (0 — ¢) + L(G, 9), (24) 
where L(@, ¢) is continuous in (6, ¢). Then, using (6), 


is 


F(@) = (29r)™' | cot (0 — ¢)f.(¢) do + 2" L( 6, &)f(o) do, 


F(6,) = N7 >’ cot 5 (6, — o.)f1(.) + wr '(2n/N) >>’ L(0, , &)f:(¢.), 


where the second integral has been approximated by a Riemann sum. Combining terms 
and noting (24), 


F(0,) = 2/N >’ (6, — 6.) "Si@.) 


=(2/r) >’ (—3s) "f(z, 2, = 8L/N. 
Nes<N 


(25) 
Equation (25) relates to (21). We conjecture that the modified quadrature formula 


T(f(a)Jenr, = (2/m) DY’ (— 8) "f(a, 


where 2, 


(26) 
= sAx, would be more suitable for Eq. (20). In fact, if one puts 


f(x) = exp (7kz), 


then one can write the summation in Eq. (26) in the form 


, , 1 + exp (—7k Az) ) 

—exp (ikz, Zia” 's k Ax) = —i exp (ikz,)9§ log - an Pn 27 
exp (ikz,) p> _ 2io sin (ok Ax) i exp (ikz,) {log i~ anton eal (27) 

Here the symbol g represents “imaginary part of” and the log function is suitably 


defined as a single-valued function. When k Az is not a multiple of z, then simplification 
of (27) and substitution into (26) gives the result 


T (f(x) leer, = @ exp (tkz,), (28) 
which is the exact result. 
More generally, consider functions of the form 


K 
f(a) = | g(k) exp (ike) dk (29) 
J -K 
where g(k) is a continuous function and K is a positive quantity less than 7/Az. The 
infinite series 


2 (r — s) exp (tkz,), 


zs, =f 2, 
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summed over odd s, is uniformly convergent in k because of Abel’s lemma. 
Hence 
»K .K 
(2/r) >’ (r — 8)" | g(k) exp (ikx,) dk = (2/r) g(k)[ >>’ (r — s)™' exp (ikax,)] dk 


K 
°K 
=1 | g(k) exp (ikx,) dk (30) 
in view of (28). However, the transform of f(x), Eq. (29) is given by 
i »K 


T(f(x)]) = 9 | (x — y)™ | g(k) exp (iky) dk dy. 
J—o@ J-K 


When x is restricted to any finite interval, the order of integration may be interchanged 
provided that one interprets the singular integral with respect to y as a “principal 


value’’. Hence 


0K 7 -) 
T(f(x)] = 3 | g(k) | (2 — y)' exp (iky) dy dk 
JK J -« 
0K 
=1 | g(k) exp (ikx) dk 
J-K 


which agrees with (30). Therefore the quadrature formula (26) is exact for functions 
of the form (29) provided that 
K Ar <x. 
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NOTE ON THE EXISTENCE AND DETERMINATION OF A VECTOR 
POTENTIAL* 
By A. F. STEVENSON (Wayne University) 
The problem of finding a solution of the equation 

curl F = f, (1) 
where f is a given vector function of position and F is to be found, is a classical one, 
and one would hardly think that there was anything fresh to be said about the matter 
to-day. In connection with some recent work, however, the writer had occasion to look 
into the matter somewhat closely, and as a result of this it appeared that the usual 
treatments’ are not satisfactory if f is only defined in a restricted region. Usually, in 





*Received June 23, 1953. 

1The following textbooks were consulted: R. Courant, Differential and integral calculus, vol. II 
(Nordeman, New York, 1944); E. Goursat, Mathematical analysis, vol. I (Ginn, New York, 1904); 
E. Picard, Traité d’analyse, vol. I (Gauthier-Villars, Paris, 1922); C. E. Weatherburn, Advanced vector 
Bell, London, 1928); H. Lass, Vector and tensor analysis (McGraw-Hill, New York, 1950). 





analysis 








1954] A. F. STEVENSON 195 
fact, either the region is not mentioned at all, or else it is assumed that it is the whole 
of space’; and the solution ordinarily given does not actually furnish an acceptable 
solution in some cases. 
It is commonly stated that the necessary and sufficient condition for (1) to possess 
a solution is 
div f = 0, (2) 
but we shall see that this may not be sufficient and that other conditions may be neces- 
sary. 
Taking Cartesian coordinates—as we shall throughout—the solution usually given 
is obtained by supposing that one of the components of F, say F,, is zero, leading to 
- 


F. = | f, dz + (2, y), Fj=- | f.dz+ (a, y), 


ve 


(3) 


where c is some constant, and, in the integrations with respect to z, x and y are held 
constant. Taking account of (2), we then find that this is a solution of (1) provided that 


the functions ¢, y satisfy 


0 0 : , , 
= - — = JAX, y, Cc); (4) 
Ox oy : 

and this equation can be satisfied, for instance, by taking y = 0, @ being then given 


by a quadrature. More generally, we could take c in (3) to be a function of x, y, or even 
take c to be different functions of x, y in the two integrals; this only leads to a modification 


of the right-hand side of (4). 

The solution (3) implies, however, that the points (zx, y, z) and (2, y, c) can be joined 
by a straight line lying entirely in the region within which f is defined. This will not be 
the case for some regions if c is a constant or any continuous function of x, y. It can 
always be achieved by choosing c¢ differently in different portions of the region, but 
this may cause the solution to be discontinuous. Similar difficulties may arise from the 
quadratures necessary to solve (4). We can, of course, always find a solution in the 
a portion of the region concerned, and it may happen that 
an analytic continuation in the whole region. But clearly 
does not provide an acceptable solution. We shall give an 


form (3) which is valid in 
this solution may possess 
cases can arise where (3) 


example below. 
We now consider the problem of solving (1) in some region R which, for definiteness, 


we take to be a simply-connected region bounded by n + 1 simple closed surfaces with 
continuous curvature S, , S,,--: ,S,, where Sp enclosed S, , --- , S, , these n surfaces 
being all exterior to each other. We may have n = 0, in which case the region is simply 
the interior of S, . We suppose that f possesses continuous first derivatives in R, and 
require a solution which also possesses continuous first derivatives in R. 

We shall now prove that the necessary and sufficient conditions for (1) to possess a 
solution are condition (2), together with the n + 1 conditions: 


| n-fdS=0, i=0,1,---,2, (5) 


Courant (Ref. 1) considers a region which is the interior of a parallelepiped, thus avoiding the diffi- 


culties mentioned here. 








~~ 
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where n denotes a unit vector along the normal drawn outward from FR at any one of 
the bounding surfaces. It follows from (2) that if any n of the conditions (5) are satisfied, 
then the remaining condition is also satisfied. If n = 0, the single condition (5) follows 
1 general, however, (5) is not a consequence of (2). That conditions (5) are 


from (2). l 


necessary is immediate, since, by Stokes’ theorem, 


S 


| n-curl FdS = 0, f= Oi, +>: ,m: 


We shall show that (2) and (5) are sufficient conditions by constructing an explicit 
solution in terms of the solutions of certain Neumann problems, using for this purpose 
a generalization of a method used in electromagnetic theory, where, however, it is always 
supposed that R is the whole of space. 

We first extend the definition of f throughout space. Inside S;(i = 1, --- , n) we 


extend it by requiring that (2) be satisfied and that n-f be continuous on S; . This can 


be accomplished by putting, inside S, 
f= VV, 
where V is a harmonic function which satisfies the boundary condition 
n-VV; =n-f on S, 


This Neumann problem possesses a solution on account of (5). We extend f outside S, 
in the same manner, and require in addition f — 0 at infinity. This can be accomplished 
exterior Neumann problem; and we note from (5) that the extended f 
since an exterior harmonic function 


an 


by solving ¢ 
will vanish at infinity to an order 1/r° at least, 
which is such that the integral of its normal derivative over a closed surface vanishes, 
must vanish at infinity to the order 1/r° at least. 

We now assert that a solution of (1), for all interior points of R, is given by 


F = curlA, A = (1/47) | f/r dr, (6) 
where r denotes the distance from a point in the domain of integration to the point at 
which A is estimated, the integral being taken throughout space. In the first place, the 
integral is convergent on account of the behavior of f at infinity. Further, div A = 0. 
For the contributions to div A from the various regions into which space is divided 
can, in the usual way, on taking account of (2), be expressed as integrals over the bounding 
surfaces. The integrals over the surfaces S , S, cancel on account of the boundary 
conditions imposed on the extended f. There remains an integral over the sphere at 
infinity arising from the region exterior to S, . But this also vanishes on account of the 
behavior of f at infinity. Hence div A = 0, and therefore 


curl F = curl cur] A = —V’A = f, 


as required. Also since f possesses continuous first derivatives in R, the same is true of F. 
If R is the unbounded region exterior to S, , --- , S, , and if f does not vanish at infinity 
to the required order, we first surround S, , --- , S, by a surface S, and proceed as before. 
We then obtain a solution valid in the portion of RF interior to S, , where S, can be taken 
arbitrarily large. To our particular solution we can, of course, add the gradient of an 


arbitrary scalar function. 
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In the usual physical applications, such as to the magnetic vector potential, the 
conditions (5) are satisfied as well as (2), and a solution exists. Occasionally, however, 
the roles of the electric and magnetic fields are interchanged, and the electric field in a 
charge-free region is expressed as the curl of a vector potential. The above discussion 
shows that this is not possible in, for instance, the region exterior to a closed surface 
which encloses a net charge. 

Two simple examples, (A) and (B), will illustrate the above considerations, the 
region R being taken in each case to be that exterior to a sphere of radius a whose center 
is the origin: 

(A) f=V(i7), r=@t+yt+2)™”. 
Condition (2) is satisfied, but (5) is not. Therefore there exists no solution. (There will, 


of course, exist solutions in other regions). 
(B) f= Vi(2’ — 7’)/r’). 


Conditions (2) and (5) are both satisfied. If we try to use the solution (3), the simplest 
choice for c is c = + ©, in which case (4) is satisfied by ¢ = y¥ = 0. Takingc = +o, 
say, (3) gives 
, d , ff) 
F,.=42 pe=-Y%Z PF, =0, 

oy Ox 


_fo¥.(e18_2 

I~@+yy\y 3r 3)’ 
provided at least one of x, y does not vanish; while if x = y = 0, z > a, we can use the 
limiting values of the derivatives of g as x, y — 0. We thus obtain a solution (with con- 
tinuous derivatives) valid in the portion of R for which z > 0. If, however, we extend 
this solution into the portion of R for which z < 0, it becomes singular on the line 


(7) 


x = y =0,2 < —a,since the function g in (7), and its derivatives, do not possess limits as 
x,y —Oforz < 0. It isreadily seen that any other choice of the limits c in (3) (for instance, 
takingc = + forz > 0Oandec = —o forz < 0) either gives a solution which is singular 


on the positive or negative z-axis, or else one which is discontinuous in R. Permuting 
x, y, 2 in (3) runs into similar difficulties. Thus there exists no solution of the form (3) 
in which one of the components of F vanishes. 

We can, however, easily obtain a solution in the form (6) by solving the required 
Neumann problem for the sphere. After some simplification we find: 


7 5 , 5 ¢ 1,5 
F, = yz/r’, F, = 2z/r’, F, = —2zy/r’. 


zr « 
This solution is independent of a and is valid everywhere except at the origin. In this 
particular case, we can obtain the solution more simply by writing 


sign Zfh.£ (:)| 
~~ 3 | % (2) ~ ay? \r/ | 


and treating the two terms of f separately: for the first term we use a solution of the form 
(3) with F, = 0, and for the second term one of the same form with F, = 0. But this is 
not, of course, the same as using (3) directly; and in more complicated cases it will be 
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difficult, and probably not possible, to express f as the sum of a finite number of terms 


for each of which we can use (3). 


STRESS FUNCTIONS OF MAXWELL AND MORERA* 
By WILHELM ORNSTEIN (Newark College of Engineering) 


Summary. A systematic process is devised for the derivation of the stress functions 
of Maxwell and Morera by the application of the theorem that a vector whose divergence 
vanishes is solenoidal. Symmetric and anti-symmetric matrices are established, and 
the elements of these matrices represent the stress functions of Maxwell and Morera 
respectively. 
*. rm 1 . . . 
Introduction. The procedure used by Maxwell to derive three stress functions 
representing six stress components at any point of an isotropic body is similar to that 
which Morera’ subsequently applied to the establishment of the corresponding stress 
functions. This procedure described also by Love* consisted in the choice of three stress 
components; then the substitution of these components into the equilibrium equation 
led to the remaining three components necessary to satisfy the equations of equilibrium. 
Later it was discovered by Sir Richard Southwell* that Saint Venant’s and Beltrami’s 
compatibility equations follow from Castigliano’s principle when the strain energy is 
expressed in terms of Maxwell’s and of Morera’s functions. 
Derivations. Neglecting the body forces, the equation of equilibrium is 
div T = 0, (1) 
where 7’ is the stress tensor. The three equations of equilibrium are obtained by cyclic 
interchange of xz, y, z in the equation 
OXx aXy , dXz 
——  —— f- <—— oe 0, (2) 
Ox oy Oz 
The equations of equilibrium (2) can be written as 
div A = 0, div B = 0, divC = 0 (3) 
where, because of the theorem that a vector whose divergence vanishes is a solenoidal 


vector, the following relations exist: 


A = curl F, B = curl G, C = curl H. (4) 
Consequently, the stress components may be written as follows: 
. oF OF, : OF OF; . OF, OF 
Xz = —? - —, Xxy =~ - =, Xz = — —- —"; 
Oy Oz Oz Ox Ox OY 
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. 0G, Gr _ 0G,  aG; a | 
ie= dys az’ ‘s= dz Ox ’ "= Ox dy ’ (5) 
7, 9Hs_ 9H, yy _2H,_ 9H: yy, _ Hs _ 9H 
sa = ay saz’ “y= dz dx’ a 0a oy 


Upon substituting (5) into the equilibrium conditions 

Xy = Yz; Xz = 22; Yz = Zy (6) 
and introducing, for abbreviation, the symbol T = F, + G,. + H; , the following equa- 
tions are obtained: 


aF 00. Mh =i oe 


Ox oy 
OF, , HG, —T) , He _ 
az oy ar vi 0, (7) 
OF, — T) , 9G, , Mi _ 
Ox oy + dz e. 


Equations (7), in the same manner as Eqs. (2), can be expressed as before: when the 
divergence of some vector is equal to zero, that vector is solenoidal. Thus the following 


expressions can be written: 


, , _ 9U; _ 9U2 _ 9U, _ Us _ 9U, _ aU, , 
i da jae Gi oe ~ ar’ Ay = “on dy ’ 
pa OTp Bhs » we ON _9V,_ 9V,, 
ities oy dz a-Tr="% ox’ ~ = dx dy ’ 8) 
| _ OW, aW, _ aW, _ aw, __ oW, aw, 
of oy dz’ G, = dz ex? #,—-T= ox oy 


To evaluate the expression I, the three terms on the diagonal of Eqs. (8) are added: 


F,+G,. + H; — 3r = —2T 


from which 


Te ee sa Us - W,) -32(%. - UD. (9) 
Upon substitution of (9) into (8) the following relations are obtained: 
py = 0s — es 5 3 Wa — Vi) — 555 (Ue — Wd — 35, (Vs — US) 
G, = 2 _ oh _ i (Ws — Va) — 55 Us — Wd) — S2(V,- Uy, (10) 
Hy = 22 — hs _ 32 Ws - Vd) -— 55 Us — W) — 95, (Ms - UD. 


Now, by use of (8) and ip the stresses in (5) can be determined as follows: 


ve Os, Ve 


ay tad ~ dye We + Ves (11) 
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aU, , Ww, -—- ’ ; 
Yy = 3 + 3 - _ (U; + W,), (12) 
Z Ox Ox 02 
TY, - FU, ae + 

Le = ae oe —_— oe a -oxts 9 is ‘ ) 
Z Ox” oy” Ox OY (Vi + Us) 3) 

‘ 10°U, 10U, Ws, , 10°W, ,10°V, ,10°W, 128°V, 

Age ee 6 2S Se eS ft oS Fe TS 5c. tS = 
2 dy Oz 2 oz Ox OY 2 Ox dz 2 dz Ox 2 dy Oz 2 O02 
cas #4. ee ee te ws 

= =-—j-—--—(V,+ U.) +=—(W, )+ — (U: - —— 
sei ye Mt to tv +20,+ Wy} ci 

: . 19°(U. + V;) 10°(W. + V3) 10(U,+ W,) ew, . 
Xy = Ye= 2 dz” + 2 Ox Oz + 2 Oy dz ~ Ox oy : (14) 
In a similar way one obtains 

" ‘ 10°(°W,+U;) , 10(W.+ V;) , 10°(V,+ U0.) 2V, " 
Xs = Zt = — 2 oy . 2 Oy dz * 2 Ox OY - Ox dz’ (15) 

i 10°(W. + Vs) 10°(U; + W,) 1a°(V, + U,) aw, 

5 = fi; = — — _ _ “ 7 — = — ° — 6 
vs Zy 2 Ox + 2 Ox OY T 2 Ox Oz Oy dz (16) 
The normal and tangential stress components expressed by Eqs. (11) to (16) are 
determined by the elements of a square matrix, which can be written as follows: 

|U, Uz, Us 
| 
V, ral (17) 


and when 
UO, =m; V2=x2, Ws = Xs 
there is a symmetric, or diagonal, matrix, whose elements form the Maxwell stress 
. v* . . . 7 
functions. With these elements the stress components are obtained from equations (11) 
to (16) as follows: 


Ox, , Ox Ox: , 9xs 


Xz = : Yy = —3 5 
oy” Oz d Oz a 
‘ dx. , o ‘ , ax 
Zz = Xe 4 Se Xy = Yr = -; Xs (18) 
Ox Oy Ox OY 
: ‘ a : : 0x» 
Ye = Zy = —-—_ Xe = Ze = = —*= 
Oy Oz Oz OX 


When the diagonal elements of the Matrix (17) are set equal to zero; i.e., when 


U,=V.=W;=0 
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and when 

W, + V; _ —V 

U; + W, = — yr 

V;, a U2 = — wz 


the matrix is anti-symmetric, and its elements represent the Morera stress-functions 
which when substituted into Eqs. (11) to (16) yield the following stress components: 


- _ Fh , _ Ove _ bs 
Xt=s ae’ 'Y™ sean? = sas’ 
ry = Yr = ~ 12 (d%: ots _ 2t3) 
aes ts = 32 (2% + oy dz /’ 
(19) 

ne eS (2% de 24s) 

is = Z¢= 2 dy \ dx oy ° dz /’ 

1 


12 (_ ads 4 Ot 4 d¥s) 
2 dz Ox oy dz /° 

Thus, the stress tensor expressed by Maxwell and Morera functions was derived by a 
direct method from diagonal and antisymmetric matrices. 


A REMARK ON INTEGRAL INVARIANTS* 
By H. D. BLOCK (University of Minnesota) 


Let the 2n variables q, , G2, *** » Qn» Pi» D2» *** » Pn be related to the 2n variables 


Q,,Q.,-:+,Q,,P,,P., ++: , P, by a canonical transformation. Let o be the unit 
square:0 S u $ 1,0 Sv S 1, and let gq; = f,(u, v), p; = gi(u, v), @ = 1,2, --- ,n), 
where f,; and g, have continuous derivatives on o. This induces the relationships Q; = F; 
(u,v), P; = G,(u, v), (@ = 1,2, --- , n). Lets; = « Sm (f;(u, v), g:(u, v)) and S; = 


iu.svee (F:(u, v), G:(u, v)), ie. the maps of ¢ on the (q,; , p;) and (Q; , P;) planes re- 


spectively. Let 


> I/ dq; dp; and > I/ dQ; dP; 
i=] ve i=l 


be denoted respectively by 


I| >> dq; dp; and I >> dQ; dP; . 
It is widely’ believed that under the conditions stated 


I > de, dp, = I > dQ, dP, . 


Cf. e.g. Goldstein Classical Mechanics, Addison Wesley, 1950, pp. 247-250. Also Corben and 
Stehle Classical Mechanics, Wiley, 1950, pp. 292-3. 
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It is the purpose of this paper, (1) to show by a simple example that this is not true, and, 
(2) to discuss the sources of the mistake and point out correct formulations. 


1. Example. Let Q; = q, , Q. = po, Pi = pi — 2p. , P2 = —2Q, —q2 . It is easily 
verified that this is a canonical transformation. Let gq, = g. = u and p, = p. = 2, so 
that Q, = u, P; = —v,Q,. =v, P, = —3u. Then s, isthe squareO Sg, $ 1,0 S p, £1; 
8 is the square 0 S gq, $ 1,0 S p. S 1; S, is the square:0 S$ Q, S$ 1, -—1 S P, £0; 
and S, is the rectangle 0 < Q, < 1, —3 S P, S 0. Then 


i > da; dp; =~ 1+ 1 = 2, 
while , 


I/ > dQ; dP; =1+3 = 4. 


Ss 


“ 


2. Sources of the mistake. “Proofs” of the incorrect statement may be made in 
several ways, two of which we now examine. First” one can show correctly that 


Og; OD; Op; Oa; ~ (dQ; oP; OP; dQ; 
(ado. _ So. 8n) 3 (20. 8B. _ 8°. 20), 


‘=1 \Ou Ov Ou ov ‘21 \Ou_ Ov Ou ov 


Now, for each 7 in (1, 2, «++ , n), 
al 


ep a1 9( - ) | 
iT dq; dp; = | | age du dv, 


and 


‘ 1 pl 9/ 2 
| a(Q, Ps du dv. 


IJ aQ, dP, = ] d(u, v) 


/0 


The error is now introduced by neglecting the absolute value signs in the integrand. 
Summing on 7 completes the ‘‘proof’’. [Since vertical lines are used to denote both de- 
terminants and absolute values and since the printing of both symbols: is widely 
used with a different meaning, it would be highly desirable if the formula for trans- 


forming of integrals were written, say, 


9 r | ee Yy) 3 
[| dx dy = | Au, >) du dv 


4 A* 
in order to avoid such mistakes. ] 
A correct formula would therefore be 


> con Gi » Pad i = ¥ son Wi Pi) I , 
> sgn Page | dq; dp; = 2 sgn age wal | dQ; dP; . 


i=1 
Ss 
Si 


A second incorrect “proof” proceeds as follows* Suppose that A is a real interval, 


2Cf. Goldstein, loc. cit. 
8This notation is used e.g. in Mood, Introduction to the Theory of Statistics, McGraw-Hill, 1950. 


4Cf. Corben and Stehle, loc. cit. 
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say 0 <A S 1. Let gq; = ¢; (A), p; = ¥; (A), @ = 1, 2, --- , n), where ¢; , y; have con- 
tinuous derivatives on A and (¢; , ¥;) traces out a simple closed curve c, as \ traverses A. 
This induces the relationships Q; = ®; (A), P; = WV; (A). Because the transformation is 
canonical, (®,; (A), Y,; (A)) also traces out a simple closed curve C; in the (Q; , P;) plane 
as \ traverses A, and 


n 


> l¢ P; dQ; - $ p. an. | = > lf W,(A) d&(A) — W(r) ae. | 


i=1 


al 


| dH(q(), pd) = 0. 


Thus, correctly, 


i=1 « 


S$ Pid = Lh rida. 


Now the mistake is made of using Stokes’ theorem £-xdy = [f,4 dx dy to complete 
Sut for this formula to hold, C must be positively oriented. Now in the 
case at hand C; and ec; have the orientations which are given to them by the mapping 


from A, and these are not necessarily all the same (examine the Example). We could 


the ‘‘proof”’. 


use Stokes’ theorem in the form 
I dq; dp; = (—1)” $ p. dq; and I dQ; dP; = (-—1)” $ P; dQ; , 
ae Si Ci 


where c, and p; are 0 or 1 according as c; or C; respectively is traversed negatively 
ositively as \ traverses A. Thus a correct formulation is 


- (—1)” I dq; dp; = 3 (—1)"' I dQ; dP; . 


Or | 


This subject, as well as other topics in mathematical physics, can be presented in a 
more unified and elegant fashion through the use of Cartan’s exterior differential calculus. 


This the writer hopes to do in a later paper. 


ON A GENERALIZATION OF SYNGE’S CRITERION FOR SUFFICIENT 
STABILITY OF PLANE PARALLEL FLOWS* 


By S. I. PAI (University of Maryland) 

In Ref. [1], Synge derived sufficient conditions for the stability of plane Couette flow 
and plane Poiseuille flow of incompressible fluid. In the present note, the conditions are 
generalized so that the problem of velocity profiles with point of inflection and with 
either finite or infinite boundary condition or both may be included. 





*Received August 10, 1953. This note was written when the author visited the Aeronautical Engineer- 


ing Department of the Pennsylvania State College during August of 1953, supported by ONR Contract 


Nonr 656(01) 
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The Orr-Sommerfeld equation [2] for stability of two-dimensional parallel flow of 
incompressible fluid may be written as follows: 
iaR(wl — w’’)¢ — iacRLo = LLo, (1) 


x is the coordinate in the direction of the flow; y, the coordinate 


where L = d’/dy’ — a’, 
perpendicular to x; ¢, the time; w = w(y), the velocity profile of steady flow; 


¢(y) exp {ia(x — ct)}, the disturbance stream function; c, the disturbance phase velocity 
which may be complex, i.e. c = c, + ic; ; a the disturbance wave number which is chosen 
positive, and R, the Reynolds number of steady flow. 

If the boundaries (finite or infinite) of the flow field are at y, and y. , the boundary 


conditions on the disturbance amplitude function are given by 


(yi) = 0, o(Y2) = 0, 


(2) 
¢'(y,:) = 0, ¢’(y2) = 0, 
where primes denote differentiation with respect to y. 
Multiplying (1) by ¢.dy and integrating for the range (y, , y2) we obtain 
iacR(I, + a7 It) = iaRQ + (Iz + 20°]; + aD), (3) 
where 
Io = | od, dy, J; = | o'd! dy, I, = | po’ dy, (4) 
Q= | [wo’d! + (a’w + w’’)do.] dy + | w'o'd, dy, 
and the subscript c indicates the conjugate complex quantity. The values of all the 
above integrals are finite for finite or infinite boundaries (see [3]). 
Taking the imaginary part of (3) and simplifying, we have 
pve — . pve . w'’ ~ 
| (w — c,)b'b. dy + at | (w — C,) + 572 |. dy = 0. (5) 
eS Jy: 2a 


If c, < 0, and w is assumed positive, so that the first integral is positive, the second 
integral must be negative somewhere in the range (7; , y2); thus we have 


ve 


Wrinin 
Ce > 9,2 -— Wrin 
74 6 | 
where w,,;, is a negative quantity. 
In the case of c, > 0, if c, > Wimax ; 
must then be positive somewhere in the range (y; , y2). Thus we have 


the first integral is negative; the second integral 


sf? 
1 u max 
Ce SS Wax T 


] 


2a” 


where w’/,. is a positive quantity. The value of c, must therefore lie in the following range: 


; ad 
Ww. Ww 
min max » 
i U mio < C, < > 2 u max * (6) 
a 


‘ 
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For plane Couette flow w’”’ = 0, we have 


Wmin < Cr < Wmax (6a) 


For plane Poiseuille flow w’’ = constant (negative quantity), Wain = 0 the sufficient 


conditions become 


Wenin 
De? < C, < Wimax ° (6b) 


Equations (6a) and (6b) have been found by Synge [1]. For flow without inflection 
point, the velocity of propagation of disturbances is either within a certain range (6a) 
or bounded on one side (6b). But for flow with inflection point, there is no bound for 


the velocity of propagation of disturbance (6). 
Equation (6), however, does not tell us anything directly about stability. For the 


discussion of stability, we take the real part of (3); we have 
—2ac,R(T; + °F?) = iaRI,(Q — Q.) + 212 + 20°T; + aI), (7) 


where 


T,.(Q — Q.) = Il 2 [ , w'd'd. ay |, | (8) 


“wi 


It is easily seen that 


| TQ — Q.) | < 2qlol: , (9) 
where g = | W’ |max - 
For any real constant a and b, we have 
M = (6 + aw’d’ + bd’’)(o. + aw’d! + bot’) dy > 0 (10) 
or 
M I;+a’ | woo! dy + BI; — a | wo, dy, 


— 2b]? — ab [ w''o'b! dy > 0. (11) 


Jy 
Hence, if a and b are any positive numbers, 
B12 > 13(2b + abp — a’q’) + It(ap — 1), 
where 
p= |w" |mex- 
Substituting (9) and (12) into (7) we have 
oc, RVI? + a7I?) < WaRgI oI, — (2b + abp — a’q? + 20°b’) 
— [i(ap — 1 + Ba’). (14) 
Therefore the stability condition c; < 0 is satisfied provided 


(b’alhg)? < 4(2b + abp — a’q’ + 2a°b’)(ap — 1 + b’a’), (15) 
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where a and b are any positive real constants satisfying 
2b + abp — a’q’ +:-2a°b’ > 0, ap — 1+ ba’ > 0. (16) 
Lessen’s criterion [3] corresponds to p = 0, g = 1 of (16). 
We may choose some values for a and 6 to get a criterion for jet flow. If we put 
a = K/ap, b = Kq’/ap’, (17) 


where K is an arbitrary constant such that K > 1 and K’q*/p* > 1, Eq. (16) is satisfied 
for all a and (15) becomes 


2 4 2 4 6 
(Rq)’ < a a + Pe a’ — _F,. 7a + (1 = _ P at 2. (18) 
©) K"q tq / Kq K’q 


Equation (18) seems to be the criterion which is suitable for jet type flow. As a tends 
to be zero, the flow will be stable for a finite Reynolds number. 
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A note on my paper 


A FORMULA FOR AN INTEGRAL OCCURRING IN THE THEORY OF 
LINEAR SERVOMECHANISMS AND CONTROL-SYSTEMS 
Quarterly of Applied Mathematics, X, 205-213 (1952) 
By HANS BUCKNER 
Replace Dz' by Dz!, in formula (36’), replace a; by a} in formula (44), replace 2Y 
by Y in formula (45). Formula (46) is correct and corresponds to corrected formula (44). 
The author wishes to thank Mr. Kenneth Geohegan and Mr. Edwinn Kinnen for 


? 


their check of formulae (36’) and (44). 


BOOK REVIEWS 


50-100 Binomial tables. By Harry G. Romig. John Wiley & Sons, Inc., New York, and 
Chapman & Hall, Ltd., London, 1953. xxvii + 172 pp. $4.00. 
The value of the individual terms and the sum of the first z terms of [P + (1 — P)}" are given for n 
in the range 50 to 100 in steps of 5 and p in the range 0.01 to 0.50 in steps of 0.01 to six decimal places 
although the last place is doubtful. The introduction, written in a clear and concise manner, describes 


the use of the tables, its application to quality control and its relation to the ratio of the incomplete 
8-function to the complete 6-function. The tabulated data are easy to read but the nature of the entries 
prevents a more uniform tabulation. 


S. L. Levy 


emer ng 
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Stress waves in solids. By H. Kolsky. At the Clarendon Press, Oxford, 1953. x + 211 pp. 
$5.00. 


This book, the thirteenth in the Monographs on the Physics and Chemistry of Materials series, is 
a welcome addition for those interested in the propagation of stress waves in solids. The first part of 
the book treats the classical theory of stress waves in elastic materials and unites many topics which 
have only appeared in the literature. It is unfortunate that the author did not insert the results for the 
reflection and refraction of a stress wave at a plane boundary between two anisotropic materials. 

The second portion of the book treats “imperfectly elastic’ materials and contains a well written 
account of the present state of affairs of internal friction in solids. Chapter VI, Experimental Investiga- 
tion of Dynamic Elastic Properties, presents the various techniques used in the measurement of stress 
waves in solids and discusses the results of the several types of measurements. The last two chapters 
of the book treat plastic and shock waves in materials and the experimental results in this field—a 
large portion contributed by the author. 

It is the reviewer’s belief that the author has accomplished the aim stated in the introduction that 


the book should be easily read by anyone trained in mathematical physics. 
S. L. Levy 


Principles of numerical analysis. By Alston 8. Householder. McGraw-Hill Book Co., 
Inc., New York, Toronto, London, 1953. x + 274 pp. $6.00. 


This is a mathematical text which is concerned with the theoretical rather than the practical aspects 
of certain computational methods. It describes and examines the principles underlying these methods 
in an effort to provide the reader with a capital of knowledge which he can use to evaluate and develop 
new techniques. The treatment throughout has high speed digital computers in mind, but much of the 
material is relevant when hand methods are considered. The subjects treated fall conveniently into two 
categories. The first concerns the solution of equations and systems of equations; the chapter headings 
being Matrices and Linear Equations, Nonlinear Equations and Systems, and The Proper Values and Vectors 
of a Matrix. The second deals with the approximate representation of functions and the chapters here are 
entitled Interpolation, More General Methods of Approximation, and Numerical Integration and Differ- 
entiation. A short introductory chapter discusses the nature and sources of errors, and the final chapter 
is a brief description of the Monte Carlo method of estimating accuracy. Each topic is developed from an 
intermediate level and every chapter concludes with an up to date, (1952), Bibliographic section. Differ- 
ential and integral equations are not treated. In the opinion of this reviewer, the value of the book would 
be greatly improved if a few worked examples were inserted in the text. 
J. FouLKES 


Discontinuous automatic control. By Irmgard Fliigge-Lotz. Princeton University Press, 
Princeton, New Jersey, 1953. viii + 168 pp. $5.00. 


In this book the author presents, in graphical form, an analysis of discontinuously operating con- 
trollers (on-off controls). By using oblique cartesian coordinates in phase space, the trajectories of the 
motion are logarithmic spirals and a graphical analysis of the “motion” is easily plotted yielding design 
criteria. Part I treats ideal systems whereas in Part II the author investigates the effects of imperfections 
in the control mechanism. In Part III the analysis for three degrees of freedom with discontinuous 
position control is treated. 

The author has a clear concise style and has included many diagrams. The reviewer noticed one small 
error in connection with the transition points which lead to end points. He recommends this book to 
engineers working with on-off controls as well as those interested in solving non-linear differential equa- 


tions by means of piecewise continuous solutions. 
SHELDON Levy 
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Statistical theory in research. By R..L. Anderson and T. A. Bancroft. McGraw-Hill Book 
Co., Inc., New York, Toronto, London, 1952. xiv + 399 pp. $7.00. 


The book is divided into two parts: Part I. Basic Statistical Theory—150 pages. Part II. Analysis 
of Experimental Models by Least Squares—227 pages. Tables include the normal, chi-square, ‘‘t’’ and 
“F”’ distributions. 

Chapter headings in Part I are in order: probability, univariate parent population distributions, 
properties of univariate distributions, bivariate and multivariate distributions and their properties, 
derived sampling distributions and orthogonal] linear functions, derived sampling distributions: normal 
parent population, introductions to point estimation and criteria of “goodness,” principles of point 
estimation: maximum likelihood, interval estimation, tests of hypotheses, special uses of chi-square. 

Two purposes of the book as stated by the authors are to provide (1) a reference book on statistical 
theory in connection with research and (2) a textbook in statistical theory. 

The reviewer prefers to endorse the reference book objective for Part I. The beginning student of 
statistics who has a good mathematical background may be a bit confused by the statistical terminology 
of the examples while the student who is weak in mathematics will encounter trouble in the brevity of 
that part of the discussion. On the other hand Part I would certainly be a good reference for a course 
primarily aimed at the material in Part II. Part II provides an excellent reference or text for problems 
of regression and analysis of variance. 

Chapter headings in Part II are in order: regression analysis general regression model with r fixed 
variates, computational methods and methods of analysis for a general regression model, curvilinear 
regression: orthogona] polynomials, least squares for experimenta] design models, the analysis of designs 
in complete blocks, the analysis of incomplete block designs, factorial experiments, the analysis of 
covariance, variance components: all random components except the mean, analysis of data with both 
random and fixed effects (mixed model), the recovery of interblock information in incomplete-block 


designs, other topics concerning components of variance: summary of needed research. 
D. R. WHITNEY 


Complex variable theory and transform calculus with technical applications. By N. W. 


McLachlan. Cambridge University Press, New York, 1953. xi + 388 pp. $10.00. 


This book, like its first edition, contains those contributions of complex variable theory which are 
most useful in treating boundary value problems by transform techniques. The use of these techniques 
is then illustrated by the treatment of a large variety of instructive physical problems, both mechanical 
and electrical. This edition has an improved presentation and contains applications of greater techno- 
logical interest than those of the first edition. 

G. F. Carrier 


Elasticity, plasticity and structure of matter. By R. Houwink. Harren Press, Washington, 
D.C., 1953. xviii + 368 pp. $7.50. 


This new edition does not differ greatly from the 1937 original. The only changes are in the chapter 
on rubber, where the earlier discussion of rubber-like elasticity has been replaced by a presentation of 
the modern statistical theory, and in the chapter on proteins, where some new figures illustrating the 
folding of linear to globular proteins have been introduced. In the latter chapter it is surprising to find 
mentioned neither the helical] structures proposed for proteins by Pauling and others, nor the role played 


by the hydrogen bond in the stabilization of the folded structures. 
STEPHEN PRAGER 























